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I . INTRODUCTION 


A major  impediment  in  the  study  of  the  dynamics  of  a particular  sub- 
merged vehicle  is  a lack  of  accurate  functional  relationship  between  the 
motion  variables  and  the  control  inputs.  In  the  uncoupled  linear  case 
such  relationships  may  consist  of  linear  transfer  functions.  Although 
their  general  form  including  the  degrees  of  the  denominator  and  numerator 
polynominals  is  known  [1],  the  coefficients  of  these  transfer  functions 
are  often  unknown.  They  must  be  determined  either  through  analytical 
formulas  involving  hydrodynamic  coefficients  or  through  identification 
algorithms  performed  upon  experimental  flight  test  data.  The  purpose 
of  this  report  is  to  present  a new  identification  method  'GRAM  Identifier' 
and  a computer  program  for  its  application  to  flight  test  data  of  sub- 
merged vehicles. 

The  method  discussed  possesses  the  following  advantages:  a)  it  is 
noniterative  and  therefore  computationally  fast,  and  b)  it  is  noise-worthy [2]- [4] . 
Only  the  single-input,  single-output  case  is  considered  in  the  present 
report.  It  will  therefore  be  assumed  that  the  flight  test  data  consist  of 
single  input  maneuvers,  each  caused  by  the  actuation  of  a single  control 
surface  while  the  remaining  control  surfaces  are  held  at  zero  deflection. 

To  aid  the  engineer,  a computer  program  entitled  GRAM  has  been  written  that 
performs  the  necessary  computations.  The  program  is  suitable  for  analysis 
of  actual  flight  test  data  as  well  as  for  a simulation  mode.  In  the  latter 
case  the  flight  trajectories  are  first  generated,  incorporating  synthetic 
disturbance  and  measurement  noise,  and  identification  is  then  performed  on 
the  simulated  trajectories.  The  simulation  mode  is  useful  when,  for  example, 
the  approximate  transfer  functions  of  the  vehicle  are  known  (from  hydrodynamic 
computations)  and  it  is  desired  to  find  efficient  maneuvers  so  as  to  develop 
a flight  test  plan. 

The  structure  of  the  report  is  as  follows.  Section  II  presents  the 
theory  of  the  GRAM  Identifier.  Section  III  gives  a user  oriented  description 
of  the  computer  program.  The  results  of  some  case  studies,  including  those 
performed  on  actual  vehicles,  are  provided  in  Section  IV.  Appendix  A includes 
the  listing  and  flow  charts  of  the  subroutines  used  by  GRAM.  Appendix  b 
provides  a brief  discussion  on  the  equivalence  of  z-domain  and  s-doraain 
transfer  functions  and  Appendix  C deals  with  the  solution  of  a key  equation. 
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II.  GRAM  IDENTIFIER 


The  identification  problem  is  formulated  with  reference  to  Fig.  1. 

The  variable  u represents  a nonzero  input  variable  — the  stern-plane 
angle,  the  rudder  angle  or  other  control  surface  deflection.  The  cor- 
responding response  y represents  one  of  the  motion  variables  — pitch  angle, 
yaw  rate  or  some  other.  In  part  (a)  of  the  figure  is  shown  the  vehicle, 
the  instrumentation  for  the  input-output  variables,  and  the  necessary 
samplers  for  digitization  of  these  signals.  Part  (b)  of  the  figure 
provides  a discrete-time  interpretation  of  the  identification  problem, 
which  is  stated  as  follows: 

Given 

i)  the  input  and  output  measurements  v(k) , x(k) , k=l,...,K, 

ii)  the  integers  n and  r in  the  model 

y(k)  + a^y(k-l)  + ....  + a^yCk-n)  = bQu(k)  + ....  + bru(k-r)  (1) 

iii)  a statistical  description  of  the  noise  processes  w(k)  and  q(k) , 
f-fnd  the  unknown  parameters  a^  and  b.  so  that  the  model  provides  the  best 
fit  (in  some  sense)  into  the  measured  data. 

Note  that  the  quantities  ai,  bj.,  y(k)  and  u(k)  are  in  fact  to  be 
estimated.  Only  x(k)  and  v(k)  are  directly  available. 


Remarks 


The  reader  familiar  with  the  principles  of  signal  sampling  may  wish  to 
skip  these  remarks. 

• Note  that  y(k)  = y(kA),  u(k)  = u(kA),  etc.,  where  A is  the  sampling 
interval. 

• In  terms  of  the  z-transform  variable  the  relationship  of  (1)  can  be 
written  as  [5] 
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(a)  Single-input,  single-output  maneuver 
y(s)/u(s)  = H(s) 


(b)  Discrete  time  identification  problem 
y(z)/u(z)  = H(z) 

Fig.  1.  S ingle- input , single-output  identification  problem. 
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y(z)  = B(z). 

u(z)  A(z) 

b b.z^+.+bz 
o + 1 r 

. . -1  , , -n 

1 + a,  z +.  • .+az 
1 n 

bo(l-Blz_1)  . . (l-6rz_1) 

d-a^'1)  ....  (l-anz_1) 


(2a) 

(2b) 


(2c) 


The  system  described  by  (1,  or  (2)  is  stable  if  and  only  if  each 

pole  a.  satisfies  the  condition  I cx . I < 1 
i 'l 


• A discrete-time  model  of  the  form  (1)  or  (2)  retains  some  infor- 
mation about  the  original  continuous  time  system.  The  faster  the 
sampling  rate  the  greater  the  information  retained.  As  a rule  of 
thumb  the  sampling  rate  should  be  about  ten  times  the  highest 
critical  frequency  of  the  vehicle  function  H(s)  = y(s)/u(s). 

When  this  condition  is  satisfied  a correspondence  between  the 
z-domain  and  s-domain  functions  may  be  achieved.  Specif ically, 
the  equivalent  continuous- time  model  becomes  (see  also  Appendix  B) 

<50  (s+qL)  • • (s+qr)  (3) 

(s+p1)(s+p2)  . . . .(s+p^) 

- ft  In  (&i)  (or  8i=e~qi^) 

- ^ In  (a±)  (or  =e  PiA ) 


The  relationship  in  (1) , or  equivalently  (2b) , may  be  written  as 

at  £ y(z)  = BT5  u(z)  (4) 

n r 

where 

- u *1 »„) 

bt  - [b0  bj  . . . br] 


H(s) 

where 
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r T ,i  -1  ~n, 

“ t1  z 2 1 

r T ,1  -rT 

— [ 1 z • • • z ] 

The  super  T denotes  the  transpose  of  a vector  or  matrix.  Equation  (4) 
should  form  the  basis  for  modeling  the  vehicel  dynamics . However , as 
discussed  later,  the  use  of  the  vector  signals  E,  y(z),  £ru(z)  leads 
to  poor  results  in  identification.  Instead  GRAM51  relies  upon  certain 
measurement  signals  shown  in  Fig.  2 generated  by  means  of  first  order 
digital  filters.  Specifically,  use  is  made  of  the  vector  signals 


y (k)  = [yo00,  y1(k),.  . . .,  yn(k)]T 
U(k)  = [un_r(k),  . . .,  un(k)]T 

consisting  of  the  measurements  at  time  instant  kA.  They  will  be  called 
output  and  input  measurement  vectors,  respectively.  Their  z-transforms 
are  denoted  as  Y(z)  and  U(z).  It  can  then  be  shown  that 


Output  measurement  sequences:  yQ(k), yn^ 

Input  measurement  sequences:  u^_r(k),  . . • »un(k) 


Fig.  2.  Measurement  filter  system 
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(5a) 

(5b) 


and  c 


are  the  coefficients  of  the  polynomial 


n"j  _o  n _i 

Z c,.z  - n (1-Q.Z  x) 

Jl-0  1-J 

The  matrix  C is  defined  in  a manner  similar  to  C^;  it  is  in  fact 
the  (r+l)x(r+l)  dimensional  top  right  corner  submacrix  of  . 


Measurement  Filter  Theorem  [2] 

If  the  signals  y(k)  and  u(k)  satisfy  (4)  for  some  parameter  vector 
A and  B,  then  the  measurement  vectors  satisfy  the  orthogonality  condition 
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where 


a = C -1A 

(6a) 

n 

B = Cr_1B 

(6b) 

Proof  The  matrices  C and  C are  upper  triangular  about  the  cross- 
diagonal, the  latter  Raving  nonzero  entries.  Hence  these  matrices  are 
nonsingular.  We  can  therefore  rewrite  (4)  as 

(C  "1A)T  C \ y(z)  = (C  “1B)C  T?  u(z) 
n n n r r r 

Substituting  6a  and  6b  and  dividing  through  by  d(z)  one  has 


T 1 „ „T  1 „ T ^ „ 

a T7~r  C Ey(z)  - B T7TT  C E u(z)  = 0 


d (z)  n ^nJ 


d(z)  r r 


Upon  substituting  (5a)  and  (5b)  this  equation  yields 

[aT  -BT](y(z)1  = 0 

LU(2)J 

The  result  sought  by  the  theorem  is  obtained  immediately  upon  taking  the 
inverse  transform.  QED 


Corollary:  Let 


and 


A-  [aT-6T]T 

f(k)  - [YT (k)  UT(k)]T 


Synthetic  parameter  vector 
Model -measurement  vector 


then 


ATf(k)  = 0 for  all  k (7) 


Measurement  vectors 

As  stated  earlier,  the  sequences  y(k)  and  u(k)  are  not  actually 
available.  Only  x(k)  and  v(k)  are,  where 

x(k)  - y(k)  + q(k) 

v(k)  * u(k)  + w(k) 

Because  the  system  of  measurement  filters  in  Fig.  2.  is  linear, 
the  following  observation  can  now  be  made 

-7- 
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Suppose  that  instead  of  processing  y(k)  the  cascade  of  filters 
processes  output  noise  q(k).  Similarly,  let  the  lower  cascade 
of  filters  process  the  noise  sequence  w(k).  And,  let  the 
resulting  measurement  sequences  be  denoted  as  q^(k)  and  w^(k), 
respectively. 

Then 


xi(k)  = yi(k)  + qi(k) 
vi(k)  = u±(k)  + wi(k) 

where  x.(k)  and  v. (k)  are  the  data-measurement  sequences  obtained 
by  processing  x(k}  and  v(k)  by  tie  two  cascades  of  measurement 
filters. 

Let 

Q(k)  = [qQ(k),.  . . . .,qn(k)] 

W(k)  = [«n_r(k),  . .,  WR(k)] 

X(k)  = [xQ(k),  . . . . ,xn(k)] 

V(k)  = [vn_r(k),  . .,  vn(k)] 

e(k)  = [QT  (k) , WT(k)]T  noise  measurement  vector 

g(k)  = [XT(k),  VT(k)]T  data  measurement  vector 

Then 

g(k)  = f(k)+  e(k) 


For  convenience,  f(k),  e(k)  and  g(k)  will  be  called  model-measurement 
vector,  noise-measurement  vector  and  data-measurement  vector,  respectively. 
To  emphasize,  the  model-measurement  vector  f(k)  is  obtained  by  passing 
the  model  sequences  y(k)  and  u(k)  through  the  measurement  filters;  the 
noise-measurement  vector  e(k)  by  passing  the  noise  sequences  q(k)  and  v(k) 
through  the  same  filters;  and  the  data-measurement  vector  by  passing  the 
data  sequences  x(k)  and  v(k)  through  the  system  of  mesurement  filters. 

Generalized  Least-Squares  Formulation  of  the  Identification  Problem 

Given  the  data-measurement  vectors  g(k),  k=l,  . . .,  K and  the 
noise-measurement  vector  covariance 
K T 

r = e E e(k)  e (k)  (E:  expected  value  operator) 

k=l 


find  the  synthetic  parameter  vector  A that  minimizes 

K T -1 

J - Z t g (k)  - f (k) ] 1 R [g(k)-f (k) ] 

k=l 


(8) 


under  the  constraint 


XT  f (k)  = 0 
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Remark 

(The  reader  may  wish  to  skip  this  in  the  first  reading) 

If  q(k)  and  w(k)  are  stationary  white  noise  processes  with  variances 

o ^ and  o 2 respectively  and  cross-correlation  coefficient  p (which  could 
q w 

of  course  be  zero)  then 


K 

R = E (K-k+1) 
k=»l 


Q(k)QT(k)0q2 

1 

2 ^ T 

! pQ(k)W  (k)0  0 
i 4 

1 

pW(k)QT(k)0  0 
(\  w 

T 

1 

i W(k)WT(k)0 2 

1 v 

1 

(10) 


Here,  [Q(k),W(k)]  = p(k)  represents  the  measurement-vector  sequences 
resulting  from  unit  pulse  (6,  = {1,0,0,...})  stimuli  at  the  measurement 
filter  input  terminals.  We  sfiall  call  p(k)  the  pulse-measurement  vectors. 


0 Note  that  R has  the  form 
7 i 

0 

W 


[R  iO  2 

11  q 

| R12P°q 

i 2 

R0iP°  ° 

L 21  q w 

! R22°w 

where  the  matrices  R.^,  * ^21  * an<*  R22  are  ^nown  (wit^out 

the  knowledge  of  0 2,  0 2,andp).  They  are  determined  entirely 

q 

by  the  known  measurement  filters.  When  either  Oq  or  0 w is  zero 
the  matrix  R becomes  known  up  to  a scalar  multiple. 


Solution  of  the  Identification  Problem 

The  solution  X and  f(k)  which  minimize  (8)  under  the  constraint  are 
obtained  by  the  la  grange  multiplier  method: 

J*  - Z | |g(k)-f(k) | | + E vk(XXf(k))  (1 

k-1  r-1  k=l 


3J» 

9f(k) 


2R_1(g(k)-f (k))  + Vk  X = 0 
(g(k)-f(k))  = \ vk  R X 
XT  (g(k)-f(k))  = y XT  R X 
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9J* 

a\T  A f(k) 

k 


(12d) 


Equations  (12c)  and  (12d)  together  yield 


XTg(k) 

\ 

2 


k T 

-f  X R X 
XTG(k) 

XTR  X 


(12e) 


Substitution  of  (12b),  (12d)  and  (12e)  gives  the  minima  with  respect  to 
f (k)  and  v^: 


J*  * I -r-  (g(k)-f(k))TX 
k=l 


= £ -f-  g (k)X 

k=l  Z 

= ^ XTg(k)gT(k)_X 

k=l  XTRX 


By  defining  the  Gram  matrix  of  the  data-measurement  vectors  g(k) 


G = Z g(k)gT(k) 
k=l 


(13) 


we  may  write 
J* 


atga 

atra 


(14) 


The  problem  now  is  to  minimize  (14)  with  respect  to  A . This  is  quite  readily 
shown  to  be  the  eigenvector  solution  of 


(G-yR)A  =0 

corresponding  to  the  smallest  eignevalue 


(15) 
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Furthermore,  it  turns  out  that  J . . = y and  f (k)  are  given  by 

minimum  1 

T 

f(k)  - g(k)  - R X 

X RX 


Remarks 

• The  actual  solution  of  the  eigenvector  problem  in  (15) 

is  contingent  upon  the  form  of  the  R matrix.  Three  different 
cases  arise  which  are  discussed  in  Appendix  C. 

• The  standard  deviations  a and  of  the  output  and  input  noise 

sequences  are  frequently  unavailable.  Under  white  noise  assump- 
tion they  can  be  estimated  approximately  as  follows.  Suppose 
that  the  useful  signal  frequencies  are  limited  to  the  frequency  band 
[0,f  ],  then  by  high-pass  filtering  one  can  estimate  the  power- 

density  of  the  noise  in  the  region  f^  to  f^  + call  this 

density  S(f^).  Suppose  the  standard  deviation  of  the  high- 

pass  filtered  signal  is  a,  then  the  standard  deviation  of  the 
original  noise  signal  may  be  estimated  as 

o2  = 2f1S(f1)  + a2  and  s^f!>  = °2/(2fz) 

• As  mentioned  earlier,  a judicious  choice  of  measurement  filters 
can  lead  to  rapid  and  successful  identification  of  the  vehicle 
transfer  function.  The  primary  consideration  in  this  choice  is 
that  the  resulting  data-measurement  sequences  y^(k)  should  be 

as  linearly  independent  as  possible  (so  that  the  cross-correlations 
between  them  are  as  small  as  possible).  For,  if  the  measurement 
sequences  were  highly  correlated,  the  matrix  G would  becomes 
ill-conditioned.  Correspondingly,  the  solution  \ would  be 
unreliable.  For  example  if  the  measurement  filters  were  chosen 
to  have  " 0 so  as  to  have  nearly  all-pass  characteristics 

(compared  to  the  vehicle's  critical  frequencies)  than  the 
resulting  measurement  sequences  have  pair-wise  correlations 
approximating  unity. 

Another  interesting  case  worth  mentioning  is  that  when 
each  measurement  filter,  instead  of  being  chosen  as  a recursive 
first  order  digital  filter,  is  replaced  by  a unit  delay.  In 
this  event  the  formulation  coincides  with  that  considered  by 
Levin  [6].  If  the  sampling  rate  for  discretizing  the  motion 
variables  is  adequately  high,  which  is  usually  true,  the  cor- 
responding sequences  y^(k)  = y(k-i)  are  highly  correlated, 
rendering  this  choice  undesirable  [7].  Poor  identification  results 
therefore  accrue. 
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• Since  the  measurement  filters  (1-Q  )/l-Q^z  ) have  low- 

pass  frequency  characteristics  witn  unit  d.c.  gain,  a con- 
venient way  to  control  the  correlation  between  the  resulting 
signals  is  to  choose  so  that  the  mean  power  of  measurement 

sequences  diminish  in  a sensible  manner.  Call  the  mean  power 
of  the  output  of  the  ith  filter  as  p^;  the  could  be  selected 

so  that 


~ „ n-i+1  ~ 
Pi  n+1  Po 


i 1,2,...,  n 


This  choice  of  measurement  filters  has  been  implemented  in 
the  computer  program  GRAM  as  is  available  to  the  engineer  on 
a select  option  basis 


• As  discussed  earlier  the  minimum  value  achieved  by  the 

criterion  function  J is  given  by  0 , the  smallest  eigen- 
value of  the  equation  posed  in  (15).  This  value  will  be 
called  the  algorithm  error.  However,  since  the  engineer 
is  interested  really  in  the  fidelity  of  output  reconstruction 
a simple  measure  of  fidelity  may  be  used.  Let  y(k)  be  the 
reconstructed  signal  obtained  by  processing  the  measured 
input  v(k)  by  the  estimated  transfer  function  (the  true  input 
u(k)  is  used  when  simulation  mode  is  used  in  the  computer 
program  GRAM).  Then  the  percent  reconstruction  error  is 
defined  as 


K K 

ESR  = 100/  Z (y (k)  - y(k))2/  Z y2(k) 
k=l  k=l 
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III.  PROGRAM  DESCRIPTION 

The  purpose  of  this  FORTRAN  program  is  to  determine  a linear  model 
from  flight-test  data  of  a submerged  vehicle.  The  method  used  is  the  'GRAM 
Identifier'  discussed  in  Section  II.  The  computer  program  is  designed  to  work 
under  three  different  modes.  When  ISIM  = 0 analysis  of  actual  flight  test 
data  is  performed;  with  ISIM  = 1 and  1STM  = 2 flight  trajectories  are  first 
simulated  and  then  identification  is  performed. 

When  ISIM  is  either  one  or  two  the  input  data  is  generated  in  subroutine 
FILLV  where  the  type  of  control  input  is  specified  by  the  parameter  INPT.  The 
flight  trajectory  is  simulated  using  the  z-domain  vehicle  transfer  function 
when  ISIM  = 1 (subroutine  RESPON)  or  using  the  vehicle  impulse  response 
when  ISIM  = 2 (subroutine  CONVOL) . Once  the  flight  trajectories  are  generated 
the  program  uses  the  facility  of  adding  synthetic  white  Gaussian  noise  (sub- 
routine CORUPT)  to  the  simulated  flight  trajectories.  Next  the  model  ident- 
ification is  performed  based  upon  the  actual  or  simulated  flight  trajectory 
through  the  'GRAM  Identifier'  method.  The  identified  model  is  used  to  recon- 
struct a flight  trajectory  which  is  compared  to  the  actual  or  simulated  flight 
trajectory  for  analysis. 

When  ISIM  = 1 it  is  possible  also  to  simulate  a feedback  system  as 
shown  in  Fig.  3 wherein  the  vehicle  transfer  function, the  input  function 
(via  option  parameter  INPT),  the  compensator  constants  and  the  gain  constant 
are  specified  to  the  program. 

The  last  of  these  is  not  read  via  data  cards;  it  is  entered  directly  in 
the  subroutine  RESPON.  Identification  is  then  performed  upon  the  vehicle 
input-output.  For  reconstruction  one  may  use  open-loop  reconstruction  or 
closed-loop  reconstruction. 


Fig.  3.  Simulated  feedback  loop;  = COMPS(I),  A = GAIN 

The  input  data  cards  on  the  subsequent  pages  give  a description  of  all 
input  variables,  and  in  so  doing  provide  an  understanding  of  the  program  use. 
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CARD  //  1 


CARD  //  2 

Variable 

Name 

(Format) 


ISIM 

(15) 


IMRESP 

(15) 


NPULSE 

(15) 
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INPUT  DATA  CARDS 

The  first  card  is  a title  card.  Columns  1 through 
80  are  available  for  an  alphanumeric  title. 

First  option  card  which  contains  eight  variables. 

Description  Columns 


Preferred 
Value  (if  any) 


N 

Order  of  system 

1-5 

- 

(15) 

MP1 

Number  of  data  points 

6-10 

(15) 

IPLT 

Plotter  option; 

11-15 

1 

(15) 

IPLT  =0  No  plots 

= 1 Plots  only  on  line  printer 

= 2 Plots  on  printer  and  CALCOMP 
plotter 

Simulation  option  16-20 

ISIM  = 0 Performs  identification 

upon  flight  test  data 

= 1 Performs  simulation  using 

z-domain  transfer  function 
coefficients 

= 2 Performs  simulation  using 
specified  impulse  response 
h(k) . 

This  variable  is  used  only  when  ISIM=2  21-25 
it  specifies  the  type  of  impulse  response 
for  the  system  being  simulated. 

IMRESP  = 0 Impulse  response  RPULSE(k) 
read  from  cards 

IMRESP  = 1 to  5 

Synthetic  impulse  responses 
generated  (see  subroutine 
C0NV0L) 

Number  of  impulse  response  points  26-30 
-14- 
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I NORM 


This  option  allows  various  matrices 
to  be  normalized  if  INORM  = 1 


31-35 


CARD  if  3 through  card  it  IX 


2 + 2* [MP1/8]  if  I SIM  = 0 


4 + 2*N 


if  IS1M=  1 


2 + [NPULSE/8]  if  ISIM=  2 

[X]  is  the  function  that  rounds  off  to  the  nearest 
integer  greater  than  or  equal  to  X. 


These  cards  contain  different  types  of  data  which 
is  determined  according  to  the  option  ISIM.  Specif- 
ically if  ISIM  = 0 the  output  input  data  is  placed  in 
these  cards.  This  information  is  placed  on  the  cards 
in  8F10.1  fields  with  the  odd  positions  containing 
the  output  data  and  even  positions  containing  the 
input  data.  If  ISIM  = l these  cards  contain  the 
z-domain  transfer  function  coefficients.  The  number 
of  coefficients  must  equal  2N+2  and  must  be  entered 
as  follows.  The  z-domain  transfer  function  is  of  the 
form  shown  below  . 


W - ^ A 

b .+b0z  +b.z  +....+b  z 
1 2 3 n+1 


l+a-z  ^+a.  z 


•+a  .nz 
n+1 


NUMCz) 

DENOM(z) 


The  coefficients  should  be  entered  each  on  a separate 
card  in  a D22.15  format  in  the  following  order. 


1,  a^  a 3 


Vl’  V "V  ~b3 'bn+l 


If  ISIM=2  and  IMRESP=0  these  cards  contain  the 
impulse  response.  The  length  of  impulse  response 
is  determined  by  the  option  NPULSE.  Eight  data 
points  are  placed  on  each  card  in  an  8F10.0  format. 


CARD  # IX  + 1 


Second  option  card  contains  sixteen  variables 


Var iable 

Name 

(Format) 


Description 


Columns 


Preferred 

Value 


Option  to  select  a specific  input 
sequence  (used  only  when  ISIM-l  or  2) 


* 1 Impulse 

=■  2 Step 

■ 3 Doublet 

- 4 Squarewave 


NCSL  TM- 204-77 


I REM 

(12) 


INPUT  = 5 Square  wave  followed  by  exponential 
= 6 Exponential 

= 7 Periodic  impulse 

= 8 Triangular  wave 

= 9 Exponential  + decaying  sinusoid 

= 10  Random  noise 

= 11-20  shifted  functions  1-10 


Option  used  to  describe  the  order  3-4 

of  the  numerator  compared  to  the 
denominator.  This  parameter  controls 
the  numerator  of  the  z-domain  model 
transfer  function.  Specifically,  it 
limits  the  numerator  degree  in  z-^  to 
N-1REM.  For  example  if  IREM  = 1,  then 
the  model  seeks  a numerator. 

NUM(z)  = b,  + b_z_1+ +b  z“(n"1) 

1 z n 


IZTS 

(12) 


QOPT 

(12) 


This  option  determines  the  type  of  5-6 
z domain  to  s domain  transformation 
that  is  performed. 


IZTS  =0  Z domain  to  S domain 

conversion  is  not  performed. 

That  is  an  equivalent  continuous 
time  system  is  not  found 

=1  An  equivalent  continuous  time 
system  is  found  (from  the 
discrete  time  transfer  function 
fl(z)  based  on  a logarithmic 
z to  s transformation) . 


=2  An  equivalent  continuous  time 

system  is  found  (from  the  discrete 
time  transfer  function  H(z)  based 
on  a pulse  delayed  z to  s trans- 
formation. 


Option  used  to  determine  measurement  7-8 
filter  pole(s).  If  QOPT  = 0 each  of 
the  measurement  filter  poles  is  set 
equal  to  the  data  value  read  as  QSAV. 

If  Q0PT=1  the  measurement  filter  poles 
are  calculated  in  subroutine  FINDQ. 
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FDBACK 

(12) 


FDREC 


ILEVIN 


IDLY 


This  option  allows  a negative  9-10 

feedback  path  to  be  added  to 
simulate  a feedback  system  f6r 
which  the  vehicle  is  the  plant. 

FDBACK  =0  No  feedback 

= 1 Feedback  is  simulated 

FDREC  is  the  variable  to  determine  11-12 
the  type  of  reconstruction  desired. 


FDREC  = 0 Open  loop  reconstruction 

= 1 Closed  loop  reconstruction 

(Note  FDREC  must  equal  zero  if  FDBACK 
equals  zero) 


This  option  is  used  when  the  LEVIN  13-14 
identification  technique  is  desired. 


ILEVIN  = 0 Gram  identification 
technique  performed. 

= 1 Levin  identification 
technique  performed. 


Delay  introduced  on  input  numerator  15-16 

XTlTM,  . “IDLY,,  . * - (N-IREM) 

NUM(z)  = z (bi+*---+b(n+l-IREM)  Z 


) 


N Order  of  model  21-25 

(15) 

MP1  Number  of  data  points  26-30 

(15) 

NPRD  Time  scale  parameter  for  input  31-35 

(15)  signal  (useful  only  when  ISIM  = 1 or  2) 

ISKIP  This  variable  determines  the  sequence  36-40 

(15)  of  points  plotted  on  the  printer. 

If  ISKIP  = 1 every  data  point  is 
plotted  and  if  ISKIP  - 5 every  fifth 
point  is  plotted, etc. 

DELTA  Sampling  interval  41-45 

(F5.0) 
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QSAV 
(F5. 0) 


NSPQ 

(F5.0) 

NSPW 

(F5.0) 


QSAV  is  the  measurement  filter  pole(s) 
used  only  if  QOPT  = 0 and  disregarded  if 
QOPT  = 1. 


46-50  0.8  to  0.95 


Noise  to  signal  power  ratio  of  the  51-55 
output  sequence  (used  only  if  ISIM  * 1 or  2). 

Noise  to  signal  power  ratio  of  the  56-60 
input  sequence  (used  only  if  ISIM  = 

1 or  2) 


CARD  # IX  + 2 


This  card  contains  the  coefficients  of 
a first  forder  compensator  in  the  forward 
path  preceding  the  vehicle.  The  general 
form  of  the  compensator  is; 


C(z) 


COMPS (2)-  COMPS (3) z"1  cs+b. 

-1  — c+a 

1 -COMPS (l)z 


The  coefficients  are  read  in  the  form 
of  3F10.1  fields.  If  no  compensation 
is  desired  a blank  card  should  be 
inserted  in  this  position. 


END  OF  FILE  CARD 

//  Card  on  IBM  360  system. 


MEMORY 


The  total  storage  required  for  the  program  is  152K  bytes,  or  approximately 
40K  words,  on  an  IBM  360/75  computer  system.  This  will,  of  course,  change  if 
the  array  dimensions  are  changed  to  meet  the  users  test  requirements.  Presently 
the  program  can  accept  up  to  one  thousand  data  points  each  for  the  input  and 
output  signals.  The  model  sought  (or  the  simulated  model  entered  when  ISIM  = 1) 
can  be  as  high  as  ninth  order.  Under  the  third  simulation  mode  (ISIM  = 2) 
the  impulse  response  can  be  of  a length  up  to  sixty  four  data  points.  It  is 
important  to  note  that  the  square  matrices  G and  Z and  the  vectors  GAMMA,  XLAMDA, 
and  COEFF  should  have  a dimension  at  least  as  large  as  (N+N+2  ) where  N is  the 
order  of  the  model. 
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OUTPUT 


The  first  line  of  output  is  the  title  which  is  followed  by  a column  of 
program  variables  as  follows: 


STARTING  SIMULATION 


SYST  = * 0°  D'j  R = 

M ♦ l = 

s :<u 

inpt  = 

9 

IR£M  = 

1 

I ZTS  = 

,'!$PQ  = 

'ISPW  = 

>j , i 'V  i(ir 

sampling 

INTERVAL 

y paramet 

ER  = r. 

QOPT  = 

1 

i PL  T = 

FDriACK= 

F0REC= 

ILLVIN= 

n 

IDLY  = 

• 1 

1 NQRM  = 

1 

COMPS  I I ) 

~ 

The  next  portion  of  the  output  is  a plot  of  the  vehicle  input  and  output. 

On  the  left  hand  side  three  columns  of  printout  list  the  serial  number  of  the 
data  point,  the  instantaneous  value  of  the  output  and  the  instantaneous  value 
of  the  input  repectively.  In  case  a feedback  loop  is  employed  (i.e.  FDBACK  = 1) 
then  the  command  input  (i.e.  the  input  to  the  feedback  system ) is  plotted 
preceding  the  vehicle  input-output  plot.  Next  the  following  title  is  printed 

GRAM  IDENTIFIER 

The  next  output  line  lists  the  values  of  measurement  filter  poles  (Q(I)).  Note 
that  N measurement  filters  are  used  where  N is  the  system  order.  All  Q(I)  are 
equal  if  QOPT  =0.  i 


Following  the  Q parameters  is  the  listing  of  the  gram  matrix  G.  This 
is  an  NPNP2  x NPNP2  matrix  where  NPNP2  = N+N+2.  The  item  printed  next  is  the 
noise  correction  matrix  Z which  is  generated  in  the  subroutine  BUILDZ.  This 
matrix  is  also  NPNP2  x NPNP2  dimensional.  The  next  line  of  printout  is  the 
Synthetic  Coefficient  Vector  XLAMDA  which  is  generated  using  the  subroutines 
S0LVE1 , 2,  and  3.  Following  the  Synthetic  Coefficient  Vector  is  the  trans- 
formation matrix  A (again  an  NPNP2  x NPNP2  matrix)  which  is  premultiplied  with 
XLAMDA  to  obtain  the  desired  parameter  vector  GAMMA.  It  is  generated  in  the 
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subroutine  BUILDA.  The  next  value  printed  is  the  estimate  of  bias  in  data 
followed  by  the  variable  NN.  At  this  point  in  the  output  the  following  are 
printed 

a)  the  z-domain  denominator,  numerator  and  poles, 

b)  the  s-domain  poles,  numerator  constants  of  a partial  fraction 
expansion  (available  only  when  IZTS  = 2 or  3),  the  denominator  and 
numerator. 

The  values  of  the  denominator  and  numerator  coefficients  for  both  the 
Z and  S-domain  are  printed  in  ascending  order  of  the  degree  of  the  term  it 
multiplies,  starting  with  the  constant  term  and  ascending  to  the  appropriate 
highest  order  term.  At  this  point  the  reconstruction  from  the  model  is 
obtained  via  RESPON.  The  subroutine  ERROR  calculates  and  prints  the  re- 
construction error 


PER  CENT  MEAN  POWER  ERROR  OF  RECONSTRUCTION  C.CJt' 

PER  CENT  OF  SQUARE  ROOT  OF  POWER  ERROR  IN  RECUSTRUCTIOM  r.0,>2 


The  last  output  is  the  plot  of  the  true  response  (when  ISIM  = 1 or  2)  or 
actual  flight  test  data  when  ISIM  = 0 and  the  reconstructed  response.  The 
same  format  is  used  as  the  previous  one  for  the  vehicle  input  output  plot. 
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PROGRAM  GRAM 


FLOW  CHART: 


START 


r — ■ 
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1 


C PROCRAM  CRAM 

01  MENS t ON  XI 1000  I , VI  10Q01 ,XOUG( ICOO ) , VORGl  100<J ) , XRECI1000I 
DIMENSION  OAT  A I 1000, 2 > ,DATA2< ICOO, 21  BUFFI  3000) 

DIMENSION  Gl?0,20)  . Z ( 20, 2C I , GAMMA  I 20 i , XL AMOAI 20 ) .COEFFI 20  ) 
DIMENSION  MPULSEI6',) 

DIMENSION  TITLEIOO) 

COMMON  NN 

REAL«B  G.Z,  GAMMA, XL AHDA . COEFF , CCMPS 

RE AL*8  DELTA,0,QSAV,DELSAV,AVGQ,AVGW,SUMV2,XSAV 

REAL  NSPQ.NSPW 

INTEGER  OOPT , FODACK , FDREC 

CQU  1 VALtNCE  I Z(l,l ), BUFF (1) ),IGI 1,11 .BUFFI  1501) ) 

EQUIVALENCE  ( OATAI 1 , I ) , X I 1 1) , I DATAU , 2) , VI 1) I 
EQUIVALENCE  I XORGU),  OAT  *7.11,1.)  I , I XRECI 1 > .DAT A2 U , 2 > ) 
EQUIVALENCE  1 NSPQ, SI GQ I , I NSPW , S I GW) 

COMMON  /COMPEN/COMPSI 1C) 

COMMON  /GKRO/IGKR 

c 

WR1TEI6, 10221 
READI5,1P21)TITLE 
WRITEI6.1021ITITLE 
HR  I TEI 6, 1023) 

MAXPI.*1000 
MAX =20 

A 320  RE ADI  5 ,1001 ) N, MPl , I PLT , I S I H, I MRESP.NPULSE, INARM 

NPNP2*N«N*2 
RDEL=0.01 

IF  I ISIM.EQ.0)REA0(5,6995) I XORGI K ) , K = 1 ,MPI ) , ( VORGI K ) , K=1 , MP1 ) 
IFIISIM.EQ.2.AN0. I MR OSP.EO.OIRE AD  15,6995)1 HPULSEC  K ),K*1,N PULSE) 
IF < ISIM.NE.1IG0  TO  6622 
C READ  DIFFERENCE  EOUATION  PARAMETERS 

RE A0I5, 701, ENO= 1234)1 COEFFI ll,l«l,NPNP2) 

CALI  FILLVl VORG.MPl.O.NPRD) 

CALL  RESPON(X,VORG,N,COEFF,XLAMOA,MPl,0) 

IF ( I PLT  .NE.  2)  GO  TO  6622 

CALL  PL0P8IHPl,l,X,MAXPL,C.0,R0EL,3, 

117HIMPULSE  RFSPONSE_, 

216HT1ME  IN  SECONDS,, BUFF) 

6622  CONTINUE 
KKKK=0 
C 

c 

4321  READ(5.l,END=123'i)INPT,IREM, I ZTS , OOPT , FDBACK, FDREC , I LEVIN , IDLY, 
IN  DEI.'  ,MP1,NPRD,  I SK I P, DELTA ,QSAV , NSPQ ,NSPW 
RE ADIS, 6995)1  COMP SII), 1*1, 31 
I F < N. EO. 10000 ) GO  TO  4320 
O'QSAV 

WRITEI6,1000IN,MPl, INPT, I RFK, I ZTS, NSPQ.NSPW, DELTA, 0, COPT, 

1IPLT, FDBACK, FOR EC, I LEV  IN, IGKR.IDLY, I NORM , 1X0MPSI l) ,1  = 1,3) 
NH1=N-1 
NP1«N*1 
NP2-N+2 
NPNPl«N*N*l 
NPNP2»N*N«-2 
RHO'O.O 
NN-N-1REM 
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IF! ISIM.EQ.OIGO  TO  23 
C 

C FILLING  the  input  ARRAY  ACCORDING  TO  OPTION  PARAMETER,  INPT 

CALL  FILLVl VORG, MPl, INPT, NPROi 
IF< FORACK.EO.OIGO  TO  6626 
00221*1, MP1 
V<n*VORG(Ii 

22  XI 11*0.0 
IFIKKKK.NE.O)  GO  TO  6616 

CALL  PtOTITIOATA  , 2, MP1, 1 ,MP1 , I SKIP , PAXPL , 1 , 1.0 » 

6616  CONTINUE 

IFIIPLT.NE.2  .OR.  KKKK.GT.OI  GO  TO  6626 
CALL  PL0PPIMPI,2,0ATA  ,MAXPL ,0.C ,ROEL , 3, 

125H  INPUT  TO  FEEDBACK  SYSTEM., 

216HTIME  IN  SECONDS., BUFF ) 

6626  CONTINUE 
C 

C GENERATING  SEQUENCE  X(KI 

I F € ISIM.E0.1ICALL  KESPCNI XORG , VORG, N ,COF FF , XL AMOA, HP  1 ,F DO AGIO 
IF( I SIM. EO. 21  CALL  CON VCL ! HPULSE , VORG, XORG, NPULSE, MPl , IMRE SP ) 

23  D0261=l,MPl 
VI  I»*VORG( I ) 

26  X(I)*XORG(I) 

IF  I (NSPQ*NSPW)*ISIM  .NE.OICALL  CORUPTl X, V, SIGQ , SIGH, MPII 
IFIKKKK.NE.O)  GO  TO  6611 
HR  I TEI 6, 100  3) 

CALL  PLOT  IT IOAT  A , 2 , MP 1, 1 ,MPl , I SK IP , K AXP L , 1 , l . 0 1 
6611  CONTINUE 

IF(IPLT.NE.2  .OR.  KKKK.GT.OI  GC  TO  6633 
CALL  PL  OP 81 MP 1 , 2, DATA  ,MA  XPL , 0. 0 ,ROEL ,3, 

127HC0RRUPTE0  INPUT  ANO  OUTPUT., 

216HTIME  IN  SECONDS., PUFF) 

6633  CONTINUE 
N-NDEN 

C START  IDENTIFICATION  FROM  INPUT  OUTPUT  DATA 
C 

CALL  GRAMII(X,V,MP1,SIGQ,SIGH,RHG,N,0ELTA,Q,Q0PT,IREP,I2TS,GAMMA, 
1 XL AMOA iG«2,MAX, I LEV  IN, IDLY, INORM) 

IFD*C 

IF ( FDBACK.GE. 1 .AND.  FDREC.EQ.l)  I FD  = F DBACK 

IF(IFD.NE.O)  CALL  F I LL  VI  VORG,  PP1 , INPT  .NPP.O  I 

CALL  ERROR! XREC, VORG, GAMMA, MPl, N, XL AMOA, XORG, I FD, IDLY) 

C 

C PLOT  RECONSTRUCTION 

HR  I TE( 6 , 66 ) 

HR  I TE! 6 , 10061 

IF!  IPLT.EQ.O)  GO  TO  6566 

CALL  PLOT1TIOATA2,2,MP1,1,MP1,I SKIP.MAXPL,) ,1.0) 

6566  CONTINUE 

99  RDEL=DELT  A 
IFIIPLT.NE.  2)  GO  TO  6666 

CALL  PLOP8(MPl,2,OATA2,MAXPl,C.O,RDEL,3, 

160MTRUE  RFSPONSE  VS  RECONSTRUCTED  RESPONSE., 

216HTIME  IN  SECUNQS., BUFF  I 
6666  CONTINUE 
C 

KKKK*KKKK«l 

100  GO  TO  6321 

1236  CALL  PICSI7<0.0.0.01 

1022  FORMAT !////// ///////// //////////////////////////////////////////) 
1021  FORMAT 180A1) 
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I 


1023 

1001 

701 

1 

1003 

1004 
66 

6995 

1000 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

0 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


FORMAT! /tlX, •*•*♦♦♦*.••♦*♦**•••*♦*♦•♦♦*•♦••*•♦•♦•*••♦*•••♦*••••, 

FORMAT (015) 

FORMAT  I 022.  19 ) 

F0RMAT(8I2,4X,U6,<.r5.CI 

FORMAT!//, SX, 'OUTPUT  !*)  ANO  INPUT  (•)',/) 

FORMAT!//, 5X, 'OUTPUT  (*)  AND  RECONSTRUCTION  (O'./l 
FORMAT!//, IX, 'TRUF  RESPONSE  VERSES  RECONSTRUCT  0 RESPONSE',//) 
F0RMATI8F10.P) 

FORMAT (1H1.50X, 'STARTING  SI MULAT I ON *,/20X,' SYSTEM  ORDER  * ',15,/, 
120X,'M  ♦ 1 « ', IS, ///,20X,' INPT  * ' i I 5 • / ,2GXi 

2 • I REM  * • , 1 5* /• 20X, • I Z TS  = • , I 5 , /// , 20X, 'NSPO  * • , F10.6 , / ,20X  , 
3'NSPW  « '.Fin. 6,///, PO*. 'SAMPLING  INTERVAL  = • , F10. 6 , /?D*  . • 3 PARAP 
4ETFR  » •,F10.6./,20X,'C0PT  * • , 1 5 , / • 20X , • I PLT  * ',I5,/,20X, 
5'FDOACK*', !5,/,20X, 'FORFC*  • , 15 , /20X, • ILEVI N*' , I 5 , < 20X , • I GKR  * ', 
615, /20X, 'IDLY  ■ ' , 1 5, //OX , ' I NORM*  • , 1 5, /20X , • CUPPS! 11  = '.3G17.10. 
7/1 
STOP 


DEFINITION  OF  PARAMETERS  UStO  IN  THE  SIMULATION  OF  A 

linear  dynamic  system 

X IS  the  corrupted  output  sequence 
V is  the  corrupted  input  sequence 

GAMMA  IS  THE  COEFFICIENT  VECTOR 

MAX  * ACTUAL  DIMENSION  SIZE  OF  2-DIM  ARRAYS  IN  THE  DIMENSION 
STATEMENT 

N = ORDER  OF  SYSTEM 

THE  MAXIMUM  VALUE  OF  N IS  MAX/2-1 

MPl  « MU,  THE  TOTAL  NUMBER  OF  SAMPLED  POINTS  IN  EACH  SEQUENCE 

SIGQ  ■ THE  STANDARO  DEVIATION  OF  THE  OUTPUT  NOISE  St QUENCE , 0! K J 
SIGH  « THE  STANDARD  DEVIATION  OF  THE  INPUT  NOISE  SEQUENCE,  W(K) 
HOWEVER  IN  THE  READ  STATEMENT  THE  DESIRED  NOISE  TO  SIGNAL  POWER  RATIO 
*******  **************** *********** 

IS  RE AO  INTO  SIGQ  AND  S1GW  FROM  WHICH  THE  TRUE  STANOARC  DEVIATIONS 
ARE  COMPUTEO  AND  STORED  BACK  INTO  SIGQ  ANO  SIGW 
RHO  * EXPECTATION!  W!KI*Q(K)  I 

DELTA  IS  THE  SAMPLING  INTERVAL 

0  IS  THE  DENOMINATOR  PARAMETER  OF  THE  KNOWN  FIRST  ORDER  0IG1TAL 
FILTERS  FOR  THE  GRAM  II  TECHNIQUE 
QSAV  IS  THEIR  CUTOFF  FREQUENCY 


IGKR-0  USE  IS  MADE  OF  THE  FIRST  ROW  OF  ADJOINT 

1 DIAGONAL  (NEGATIVE  ENTRIES  SET  TO  ZER01 

2 ABSOLUTE  VALUE  OF  0 I AGONAL 


I SI M*0  READ  EXPERIMENTALLY  MEASURED  INPUT  OUTPUT  DATA 

1 REAO  COEFFICIENTS  OF  H(Z>,  THEN  SIMULATE  INPUT-0U1PUT 

2 READ  IMPULSE  RESPONSE  HPULSEm,  K«l, .. . .NPULSE  ANO  SIMULATE 


IREM  - DEGREE  OF  DENOMINATOR  MINUS  DEGREE  OF  NUMERATOR 


• "25" 


i 


i 


oor»oooooooooooooooooooooor»oooooooooooooooooooooor»r*oooo 
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IDLY  « DELAY  INTRODUCED  ON  INPUT  NUMERATOR;  ZETA=>1/Z 

ZETA**IDIY  ( HUM..  «BIN*l-IREM)*ZCTA**(N-IREM)  » 

IZTS  * 0 PRINTING  OF  DISCRETE  TIKE  TRANSFER  FUNCTION  ONLY  AND 
THE  POLES  OF  THE  Z DOMAIN 

» 1 IF  LOGARITHMIC  TRANSFORMATION  IS  TO  AE  CALCULATED 
* 2 II  DELAYED  PULSE  INVARIANT  TRANSTORMAT ION  IS  TO  HE  CALCULATED 

INPT=l  IMPULSE,  2:  STEP,  3S  COUPLET  (DURATION  NPUL) 

A:  SOWAVE  (PERIOD  NPUL I . S:  SQ-EXP,  6:  EXP,  7:  PRO  IMPL 

8:  TRI  WAVA,  9:  EXP40SC,  10:  RANDOM 

QOPT  - 0 IF  Oil )=OSA V 

1 IE  QIII  GENERATED  IN  QFIND 

FOR  AN  UNSTABLE  SYSTEM  FEEDBACK  MAY  BE  PROVICFD  PY  SETTING 
FEEDBACKS.  A COMPENSATOR  IN  THE  FORUM. 0 PATH  IS 
PROVIDED  ON  A NONOPTIONAL  BASIS,  EXCEPT  THAT 
WHEN  THE  COMPENSATOR  CARD  HAS  ZERO  ENTRIES  THE  PROGRAM 
AUTOMATICALLY  SETS  THE  COMPENSATOR  TO  C 1 2 1 = 1-0 
FEEDBACK  GAIN  IS  MANUALLY  ENTERED  IN  SUBROUTINE  RESPON 
AS  THE  VARIABLE  "GAIN" 

• *««»***THE  COMP S ( I ) COEFF ICI ENTS 

OF  cm=  <C0HPS(2I-CCMPS(  ■)}//.  ) / (1  -COMPSID/Z)  MUST 
MUST  BE  READ,  A BIANK  CAPO  MAY  BE  PROVICEO  IF 
IF  NO  COMPENSATOR  IS  CESIRED 
FDREC  =0  IF  OPEN  LOOP  RECONSTRUCTION  IS  DESIRED 
MUST  BE  ZERO  IF  FDBACK  IS  ZERO 
1 IF  CLOSCO  LOOP  RECONSTRUCTION  OESIREO 


POLES  OF  THE  CONTINUOUS  DOMAIN  MUST  BE  DISTINCT  AND  NON-ZERO 
FOR  TRANSFORMATION  TO  BE  VALIC 

IT  IS  IMPORTANT  TO  NOTE  THAT  THE  VALUES  OF  THE  CONTINUOUS 
SIGNALS  SAMPLED  AT  T I ME  - I K-l ) +DEL  T A ARE  STORED  IN  THE  KTH 
SEQUENCE  POSITION  OF  THE  ARRAYS 

DATA  DECK  CONSISTS  OF  A READ  OPTION  CARD  IN.IPLT), 

THE  Z-DOMAIN  COEFFICIENTS  OF  THE  ORIGINAL  TRANSTFR  FUNCTION, 
AND  A PARAMETER  CARPI N , MP 1 , 1 NPT , I GRM , 1 RE M, l STZ , S I GQ , S 1 GW , 
D£LTA,QSAV,Q0PT,NPRC,F08ACK,FCREC). 

PROGRAM  READS  SEVERAL  PARAMETER  CARDS.  LAST  PARAMETER  CARO 
MUST  HAVE  1 IN  COL  1 TO  READ  ANOTHER  TRANSFER  FUNCTION  ANO 
PARAMETER  CARD  SET. 

I PLT=0  NO  PLOTS 

IPLT=l  PLOTS  ONLY  WITH  PRINTER 

I PLT*2  PLOTS  ON  CALCOMP  AS  WELL  AS  PRINTER 


ENO 
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IV.  APPLICATION  EXAMPLES 

Example  1 


For  a six-man  submersible  the  transfer  functions  describing  its 
dynamics  were  obtained  from  the  hydrodynamic  coefficients  via  a 
computer  program  RGEORGF. . The  pitch  vs.  stern-plane  dynamics  will  be 
used  here  to  demonstrate  the  application  of  GRAM  Identifier  (ILEVIN  = 0 
in  the  program).  Specifically  the  transfer  function  relating  these 
variables  is 

„ . . -0.3A320s2  -0.17384s  -0.008631 

H4^  4 3 2 

s + 1.47989s  + 0.11833s  + 0.02048s  + 0.00102 

-0.08348  + j0. 4836b  -0. 08348-1 0.48366 

s + 0.00919  + jO. 11378  ' s + 0. 00919-J0. 11378 

0.16696 0.00002 

s + 1.4057  s + 0.05579 


For  all  practical  purposes  this  is  seen  to  be  equivalent  to 


H4(s)=  H3(s) 


-0.34320s  -0.15469 

s3  + 1.58950s2  + 0.26054s  + 0.00306 

-0.08349  + j0. 48367  -0.08349-10.48367  0.16696 

s+0 .00919+10. 11378  s+0.00919-j0. 11378  s+1.4057 


In  this  third  order  function  the  energy  of  the  complex  pole  pair  is  26.9 
while  the  energy  associated  with  the  real  pole  is  0.01;  the  latter  thus 
represents  only  0,037%  of  the  energy  at  the  dominant  complex  pole  pair. 
Therefore,  unless  the  input  is  such  that  its  spectral  content  is  rich  in 
radian  frequencies  around  1.4,  this  mode  will,  be  extremely  feeble.  We  will 
in  fact  call  this  mode  a micromode  [3].  When  this  micromode  is  not 
appropriately  excited  the  vehicle  transfer  function  may  be  approximated  as 

H (s)  = H„(s)  = -5 

* s + 0.1838s  + 0.01303 

-0.08349+10.48367  -0.08349-10,48367 

s+0. 00919+j0. 11378  s+0. 00919-j0. 11378 


Before  describing  the  various  experiments  conducted,  the  z-domain  description 
of  the  pitch  transfer  function  is  first  provided.  The  transfer  function 
H^(s)  was  transformed  by  the  Leading-Edge-Pulse  Equivalence  method  (Appendix 
B).  A sampling  interval  A = 0.5  second  was  used  yielding 

H (z)  (10~4)(0.20744z~1-0.19065z~2-0.15179z.~3+0.13714z~4) 

4 l-3.45527z_1+4. 38952z-2-2.41135z-3+0.47714z-4 
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In  all  of  the  five  experiments  performed  250  seconds  of  simulated  data 
(MP1  = 500)  were  used.  This  represents  approximately  2 1/2  time  constants 
of  the  dominant  mode. 

Experiment  1 

The  function  H^(z)  is  employed  to  simulate  the  discrete-time  trajectory 
0(k)  to  a given  stern-plane  input  6s(k)  (INPT  = 9,  NPUL  = 100).  As  stated 
earlier  G(s)/6  (s)  is  effectively  a third  order  function,  however,  a fourth 
order  identification  was  first  performed  without  masking  the  data  with  noise 
(NSPW  = NSPQ  = 0.0).  The  identification  yielded  the  following  H^s)  with  the 
option  variables  chosen  as  ISIM  = 1,  IREM  = 1,  IDLY  = 1,  QOPT  = 1 and  IZTS  * 2: 


H4(s)  = 


-0.34323s  -0.17 394s-0. 0086 7 

s4+l . 48029s2+0 . 11867s2+0 . 02049s+0 . 00103 


ESR  = 0.000 


(see  page  12) 


Clearly  this  identification  is  good  since  the  model  found  is  almost  identical 
to  the  given  H^(s).  A comparison  of  the  poles  of  the  given  H^(s)  and  the 
identified  poles  is  shown  below. 


H^(s)  Poles 

-0.00919  + jO. 11378 
-0.00919  - jO. 11378 
-1.4057 
-0.05579 


Identified  Poles 

-0.00919  + jO. 11378 
-0.00919  - jO. 11378 
-1.4058 
-0.05603 


Although  the  vehicle  transfer  function  was  identified  perfectly,  any 
quick  conclusions  as  to  the  effectiveness  of  the  identification  method  are 
misleading.  Because  even  the  slightest  amount  of  noise  on  the  data  will 

mask  the  micro-micro-mode  making  its  identification  impossible. 

Therefore  the  remaining  experiments  will  pertain  to  third  order  identification 
except  the  last  one.  The  latter  is  a second  order  run  demonstrating  the 
detection  of  the  dominent  pole  pair. 

Experiment  2 

This  run  is  identical  to  the  previous  one  except  that  a third  order  model 
was  sought  (N  “ 3) . The  transfer  function  found  was 

-0 . 00435s2-0 . 33603s-0 . 14986 


H3(s)  = 


s3+1.38026s2+0.03806s+0. 01774 
ESR  = 0.447 
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The  identified  poles  are  compared  to  the  true  poles  below 


H^(s)  Poles 

■0.00919  + jO. 11378 
■0.00919  - jO. 11378 
•1.4057 


Identified  Poles 

-0.00919  + jO. 11378 
-0.00919  - jO. 11378 
-1.36187 


Experiment  3 

Ten  percent  rms  noise  was  added  to  both  the  input  and  output  data 
(NSPW  = NSPQ  = (0.1)2).  A third  order  identification  was  performed  using 
an  input  generated  via  option  INPT  = 9,  NPUL  = 100.  The  test  failed  due 
to  the  extremely  poor  spectral  content  of  the  input.  Specifically  the  input 
signal  did  not  have  sufficient  energy  at  radian  frequencies  around  1.4, 
consequently  it  did  not  properly  excite  the  pole  at  that  location.  A 
different  input  was  therefore  used  (INPT  = 5,  NPUL  = 10,  QOPT  = 1).  The 
input  and  output  signals  are  shown  in  Fig.  4a.  The  corresponding  results 
yielded  by  GRAM  are  as  follows: 

~ . . -0.50521s2-2,75286s-2.5146 

^3^s^  3 2 

sJ+l. 76398s  +0.04714s+0. 02249 

ESR  =»  5.275 


A comparison  of  the  identified  poles 

Hj(s)  Poles 

-0.00919  + jO. 11378 
-0.00919  - jO. 11378 
-1.4057 


to  the  true  poles  is  given  below. 

Identified  Poles 

-0.00982  + jO. 11312 
-0.00982  - jO. 11312 
-1.74435 


This  input  signal  contained  sufficient  energy  around  the  radian 
frequency  1.4  to  excite  the  micro-mode  just  enough  for  identification 
purposes.  Figure  4b  shows  the  reconstructed  output  comparing  it  to  the 
true  output. 

Experiment  4 

This  experiment  demonstrates  the  importance  of  the  choice  of  the 
measurement  filter  pole(s) . The  measurement  filter  pole  QSAV  was  varied 
(QOPT  * 0 to  disable  automatic  filter  pole  selection)  and  its  effect  studied 
on  the  identification  algorithm.  Each  of  the  runs  performed  uses  the  options 
INPT  - 5,  NPUL  - 10,  IREM  - 1 and  IDLY  = 0 seeking  a Lhird  order  model. 

Ten  percent  noise  was  added  of  the  same  manner  as  in  Experiment  3. 

The  following  results  were  obtained: 
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CORRUPT  f.D  INPUT  UNO  OUTPUT 

I I I I I I I I I I I I I I I I I I 1 1 L-L  -LJ  Ll±-1  L_Li.  il-l-Ll  -j.-L  1_L  l.X-Ll-l-l-l-t_l-L 


a ft 
2 ‘.  .0 
l 25 
U 

-l  20 
-2  50 
-3  ft 
-0  OQ 
-G  20 
-7  50 
-6.  75 
-10  0 
-ll  J 


<1 1<  ;'i  ,'i  1 1 ii  'i  <■;  n r,  r;,  i ;j  T>  fT  s' 

ii1,!!  1!  life  ! Itili;!  !!i pi 


f > v- 


H rrrr  mT-iTi 


(a) 


TRUC  RfSPONSi:  VS  KX  CON  VIRUS  IT!)  RFSPONSf 

3 75 
2.!>0 
1.25 
0 

-1.25 
-2.50 
-3  75 
-5.00 
-5.25 
-7  SO 
-ft.  ft 
-10  0 
-11.  B 


(b) 


Figure  4.  Experiment  3 - Third  Order  pitch 
vs  Stern-Plane  Identification 

-30" 


i* 


NCSL  TM-204-77 


QSAV  (Percent 

0.70 

0.80 

0.85 

0.90 

0.95 

0.98 


ESR 

Error  to  Signal  ratio  RMS) 
00 

33.986 

14.374 

3.947 

4.556 

5.696 


This  example  should  make  the  user  aware  of  the  flexibility  provided 
to  the  test  engineer  by  the  measurement  filter  pole(s).  The  value  of  the 
measurement  filter  pole  should  be  such  that  each  successive  measurement 
filter  attenuates  the  input  signal  by  a reasonable  fraction;  in  particular 
the  output  of  the  last  measurement  filter  should  not  be  an  order  of  magnitude 
lower  in  power  than  the  input  signal  to  the  first  measurement  filter.  More 
information  on  the  choice  of  measurement  filters  is  available  in  reference  [2]. 

Experiment  5 

This  final  experiment  demonstrates  the  detection  of  the  dominant  pole 
pair.  The  input  (INPT  = 6,  NPUL  = 200)  and  output  signals  were  masked  with 
10%  rms  noise  (NSPW  = NSPQ  = (0.1)^)  and  the  following  option  parameters 
were  used:  QSAV  = 0.95  IREM  = 1,  IDLY  = 0.  The  input  used  (INPT  = 6)  is 
an  exponentially  decaying  function  whose  time  constant,  and  therefore  the 
cutoff  of  the  power  spectrum  curve,  is  controlled  by  the  value  of  NPUL 
(Time  constant  = NPUL  * DELTA).  Therefore  to  identify  the  slow  poles 
(-0.00919  t jO. 11378)  an  input  rich  in  radian  frequencies  around  0.114  is 
desired.  A vlue  of  NPUL  = 200  produces  an  adequate  power  spectrum.  The 
input  and  output  signals  are  shown  in  Fig.  5a.  The  computer  program 
identified  the  second  order  mode  as  follows 

H (s)  = 0.02308S-H). 09245 

2 s2+0.01904s+0. 01318 

ESR  = 2.641 

A comparison  of  the  identified  poles  and  the  true  poles  is  given  below. 

^(s)  Poles  Identified  Poles 

-0 . 00919+J0 . 11378  -0 . 00952+j0 . 11439 

-0.00919-j0. 11378  -0.00952-j0. 11439 

Figure  5b  shows  the  reconstructed  output  comparing  it  to  the  true  output. 

The  preceding  experiments  should  help  the  user  to  better  understand  the 
significance  of  the  option  variables  (originally  defined  on  pages  14-18) . 
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CORRUPTED  INPUT  FIND  0UTPU1 


O.  g 4CiD  g W8  J.S'J  2.Q0  2.  “9  2.ft>  3.0  3.9'J  4.«9 


(b) 

Figure  5.  Experiment  5 - Detection  of  Dominant  Complex  Poles 
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Example  2 


Flight  tests  on  a towed-sonar  vehicle  (SMS2619 one-t'nird  scale  model)  , 

designed  by  the  Naval  Coastal  Systems  Laboratory,  were  conducted  at  the 
Naval  Ship  Research  and  Development  Center,  Carderock,  Md.  The  test  data  were 
recorded  on  magnetic  tapes  at  a sampling  rate  of  30Hz.  In  this  example  the 
results  of  two  experiments  are  presented,  one  pertaining  to  pitch  vs.  stern- 
plane  deflection  and  the  other  pertaining  to  roll  vs.  rudder  deflection.  In 
both  cases  the  data  were  preprocessed  by  a 211*.  digital  filter  to  substantially 
remove  an  undesirable  3.3Hz  oscillation,  suspected  to  be  caused,  by  an  artifact 
in  instrumentation.  The  data  were  then  sampled  at  5Hz  (Nyquist  frequency  = 2.5Hz 
5 15.71  radian/sec)  for  use  by  the  identification  program. 

For  the  pitch  vs.  stern -plane  data  a second  order  (N=2)  identif ication  was 
performed  which  yielded  the  model 

u \ - 0.784(s  + 0.028) 

H2  ^(s)/&s(s)  (s  + l.  21)  (s  + 0.03) 

The  reconstructed  output  is  shown  in  Fig.  6 together  with  the  output  data 
used  for  identification.  It  is  worth  mentioning  that  higher  order  models  were 
also  attempted,  however,  the  additional  poles  found  had  insignificant  energies 
associated  with  them. 

For  the  roll  vs.  rudder  data  the  results  of  a sixth  order  (N  = 6)  identifi- 
cation are  presented.  The  model  found  is 

H (s)  = -°-09682s5+1-01434s4+0-19137s3f0-04870s2+0-00132s+J- 00013 

s6+O.67744s5+-O.6497ls4-t-0.12743s3+0. 02647s2  K>.00078s+0. 00007 


_ -0,35521  i jO. 79885  -0.00106  ± 10.00045  -0.01183  ± 0.01648 

s+0. 23747  ± 10.66789  s+0. 01006  ± 10.05316  s+0. 09119  ± jD.190L5 

The  reconstructed  output  is  shown  in  Fig.  7 together  with  the  output  data. 
However,  the  energies  associated  with  the  last  two  pairs  of  poles  are  quite 
small  so  that  a second  order  identification  would  therefore  seem  desirable. 

In  our  analysis  of  the  SMS  vehicle  we  also  conducted  multiinput-multioutput 
identification,  and  such  analysis  showed  that  the  lateral  dynamics  is  essentially 
governed  by  two  pole-pairs, one  pole  pair  reasonably  observable  in  the  roll  data 
and  the  other  in  yaw-rate  data.  The  multiinput,  multioutput  version  of  GRAM 
identifier  is  discussed  in  [9]. 
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sms  wiia  run*  'io.  piipi 

7 5c  .(-Ll.Li_Lj.UI  l1_LL.L1.Li  I LlLi  LI  1 l-LLl-l  l.l  l.LL.I-l  l.  LlLi.I  _L.l  1.LI  I 


0 13.0  JO . I *5.1  OH  ? 7b  T 35. 2 I M.  170  135. 


Fig.  6 Identification  of  longitudin.nl  dynamics  of  SMS2619 

vehicle;  reconstruction  of  pitch  from  identified  model. 


SMS  DfiTfl  RUN*  co]  . ROLL 


C.  1 V?  in  * W 0 'jA.7  70  3 r..l  9T>  » 13  l?b 


Fig.  7.  Identification  of  lateral  dynamic*  of  SMS2G19  vehicle; 
reconstruction  of  toll  from  Identified  model. 
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APPENDIX  A 


SUBROUTINES  USED  IN  PROGRAM 


The  following  subroutines  used  by  the  GRAM  Identifier  program  are  described 
in  the  appendix. 

BUILDA 

BUILDZ 

CONVOL 

CORUPT 

ERROR 

FILLV 

FINDQ 

GRAMII 

IZTOS 

POLCON 

PRCVEC 

PRMAT 

PRVEC 

RESPON 

SOLVE1 

SOLVE2 

SOLVE3 

ZTOS 

The  subroutines  DOBINV,  PLfiP8  and  POLRT  are  not  detailed.  Their  function  is 
indicated  below  and  they  can  be  substituted  by  standard  routines  from  a scientific 
package. 

DOBINV 

PLOP  8 

POLRT 


Inversion  of  a square  matrix 
X-Y  plotter  (CALCOMP)  routine 

The  subroutines  called  by  PL0P8  are  also  not  discussed  here 
Computes  the  real  and  complex  roots  of  a real  polynomial 
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SUBROUTINE: 

PURPOSE: 

EQUATIONS : 

FLOW  CHART: 


BUILD  A 

CALCULATES  THE  TRANSFORMATION  MATRIX  TO  CONVERT  SYNTHETIC 
PARAMETER  VECTOR  TO  THE  DESIRED  PARAMETER  VECTOR 


Aij  Ai,j+-1 


Q1  Ai-l,l+l 
Pj 


th 


after  1st  row  and  (n+1)  “ column 

A’JTA  ‘J 

"A^s  are  generated  by  other  means  with  i = 2,3,....,  n+1; 

3 ~a ,n-l ,...., 1 
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SUBROUTINE:  BUILD  A 

DESCRIPTION:  The  "A"  Matrix  is  formed  in  the  following  manner 


where  each  quadrant  is  an  <h+i  x (n+1) square  matrix. 

This  matrix  embodies  the  relationship  between  the  synthetic 
coefficient  vector  XLAMDA  derived  from  the  GRAM  matrix  and  the 
true  coefficient  vector  of  the  systems  transfer  function, 
GAMMA.  This  relationship  is  dependent  upon  the  values  used 
for  the  measurement  filters.  Multiplying  XLAMDA  by  a matrix 
of  order  2 (n+1)  with  r in  the  upper  left  and  lower  right 
quadrants,  yields  GAMMA. 


PROGRAM  VARIABLES: 


A "A"  MATRIX 

DEL  MEASUREMENT  FILTER  NUMERATOR 

MAX  MAXIMUM  ROWS  PERMISSIBLE 

N ORDER  OF  SYSTEM 

Q MEASUREMENT  FILTER  POLE(S) 
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SUBROUTINE  BUILOAI A.C.OEL .N.MAXI 
RE  AL*8  a( MAX, 1) ,Otl) ,OEt ( 1) .PROD 
NP1=N*. 

NPNIP2  = N + N*2 
A( l.NPi )*1.3000 
PR3D=1 . ODOL 
D0312K=1.N 
I«MP1-K 

PR0D=PR3D/DEL ( I ) 

All,n*PROD 

312  AIK+l.NPl  )=0.0DC'0 
033131=2. NP1 
D3313K=1,N 
J*MPl-< 

313  All  .J)  = ( At  J.J*1  )-0«  J)*Al  I-l.J+l  J I/DEHJ) 
D331AI *1 » NP1 

D03 1 A J = 1 » NP1 
A* I « J+NP1 ) =0. OOOC 
A< I ♦NPl, J 1=0.0000 
31A  All  AfiPl . J+NPT.  )*A(I»J) 

MR  I TEI 6 .1005) 

10J5  P3RHAT  C1X»* A- MATRIX*) 

CALL  PRMATI A.NPNP2.NPNP2.MAXI 
RETURN 
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SUBROUTINE:  BUILD  2 


PURPOSE:  CALCULATE  NOISE  CORRECTION  MATRIX  "Z" 


EQUATION : MP1 

2 R1(k)Rj(k) 


FLOW  CHART: 


R(l)*1.0 
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SUBROUTINE:  BUILD  Z 


DESCRIPTION:  This  subroutine  calculates  the  contribution  of  the  noise  to 
the  gram  matrix  resulting  from  the  measurement  filter  output. 
The  total  Z matrix  is  formed  in  four  sections.  The  First 
section  (zO-D)  is  generated  through  use  of  the  GAMMA  matrix 
The  second,  third  and  fourth  sections  deal  in  optimizing  the 
estimate  of  the  noise  correction  matrix. 

PROGRAM  VARIABLES: 


DEL  MEASUREMENT  FILTER  NUMERATOR 

ILEVIN  VALUE  IS  EITHER  0 OR  1. 

0 GRAM  TECHNIQUE  IS  PERFORMED 

1 LEVEN  TECHNIQUE  IS  PERFORMED 

MAX  DIMENSION  SIZE 

MP1  NO.  OF  DATA  POINTS 

N ORDER  OF  SYSTEM 
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SUBROUTINE: 


BUILD  Z 

Q 

R 

RHO 

SIGQ 

SIGW 

Z 

ZP 


MEASUREMENT  FILTER  POLE 
COEFFICIENT  VECTOR 
EXPECTATION  OF  (W(K)*Q(K)) 

STANDARD  DEVIATION  OF  OUTPUT  NOISE  SEQUENCE 
STANDARD  DEVIATION  OF  INPUT  NOISE  SEQUENCE 
NOISE  CORRECTION  MATRIX 
WORKING  ARRAY 
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SUBROUTINE  BH I LDZ ( Z , Z P, R . N, MP 1, S I GW , SI  GO ,RHO, OE L . 0 . MAX, I Lc VIN ) 

: SUBROUT  I Nc  FOR  CALCULATING  THE  NOISE  CORRECTION  MATRIX,  Z, 

Z FOR  GRAMI  AND  GRAMI I 

DIMENSION  Z(M4X,1 ) ,ZP! MAX, 1 1 ,R(1 ) .0(1). DELI  1 ) 

DOUBLE  PRECISION  Z , ZP , R, OEL , G, DCON 

NP 1 *N* I 

NPNP2=N*N*2 

R 1 1 ) = l , PDCC 

DOlZI'l.N 

RI  I ♦!  )=R<  I )*r»EL  1 1 » 

IF( I LEV  IN. EO. 1)R( 1*1) =0.0000 
12  CONTINUE 

D01I=1,NP1 
001 J=1,NP1 
ZP(  I,  J I = 0.0D',0 

1 Z ( I , J ) *3. ?DC" 

D02K=1,MP1 
003  1=1, NPI 
D03  J=1,I 

ZP!  I. J)  = ZP!  I, J)+RU>*R!J> 

3 Z( I , J)=Z ( I , J)+ZP( I , J) 

IF!ILEVIN.EQ.0)R!1)  = 0.  SDQ'.- 

OOA I =1 ,N 

IF!  UEV'IN.EQ.l)R!I+l)=Rm 

IFI  ILEVIN.EQ.DGO  TO  A 

RII  + l)  = RI 1 + 11*0(1 l+R( I )*DEL ( I ) 

A CONTINUE 

2 CONTINUE 
D0A0I=1,N 
IP1*I+1 
D0A0J=IP1,NP1 

ao  zii,j)=zu,n 

SI GQ2=S I GO*SI GO 
SIGW2=SIGW*STGW 
0051=1, NPI 
D05J=1,NP1 

5 ZINP1*I,NP1+J)=ZII,J)*SIGW2 
IFIRH0)6,7,6 

7 0091=1. NPI 
DOB J=1 , NPI 

ZI I+NP1,J)=0.0D00 

8 Zl  I ,NP1  + J)=0.3DC!) 

GO  TO  9 

6 00111=1, NPI 
031 1 J* 1 , NPI 

Z(I,NP1*J)=Z(I,J)»RH0*SI GW*SIGQ 
11  ZII*NP1,JI=Z(I,NP1*J) 

9 03101=1, NPI 
DOiOJ-l.NPl 

10  ZII,J)=Z(I,JI*SIG02 
MR  I TE( 6 ,1000) 

1000  FORMAT! IX, ‘NOISE  CORRECTION  MATRIX,  ZM 
CALL  PAMATIZ,NPNP2,NPNP2»MAX! 


RETURN 

ENO 
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SUBROUTINE: 

PURPOSE: 

EQUATION: 


NCSL  TM-204-77 

CONVOL 

PERFORM  CONVOLUTION  OF  HPULSE  AND  VORC 
NPULSE 

XORG  - E VORG,  , , HPULSE 

, IctI  “TO  m 

m-1 
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SUBROUTINE:  CONVOL 

DESCRIPTION:  This  subroutine  determines  the  convolution  (XORG(K))  of  VORG(I) 
and  HPULSE(J) . NPULSE  is  the  length  of  the  impulse  response 
and  the  variable  IMRESP  is  the  option  used  to  specify  the  type 
of  impulse  (HPULSE)  response  desired.  Specifically  when  IMRESP 
equals  zero  the  impulse  response  is  entered  as  data. 

PROGRAM  VARIABLES : HPULSE 

IMRESP 

MP1 

NPULSE 
VORG 
XORG 


IMPULSE  RESPONSE 

OPTION  TO  DESIGNATE  TYPE  OF  HPULSE  TO 

BE  GENERATED 

NUMBER  OF  DATA  POINTS 

NUMBER  OF  DATA  POINTS  OF  IMPULSE  RESPONSE 
CORRUPTED  INPUT  SEQUENCE 
CORRUPTED  OUTPUT  SEQUENCE 
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S J8RQUT l NE  C'TNVOLIHPULSE.  VORGrXOAG. NPULSE, PP1  , IMRESP ) 

C PERFORMS  CONVOLUTION  OF  HPULSE  AND  VORG 

C X3RG<*)=  HPULSr (1 >*V3RG(K  )♦  ....  ♦HPUlSc ( KK ) «VORG( K-KK* 1 » , 

C WHE° E KK=NPULSE  OR  K hHICHEVrft  SMALLER 

DIMENSION  HPtILSct  1 ) . VORGI  1)  , XORGI  L) 

I F ( IMRESP.NE.O) GO  TO  2: 

2 03  5 K«X,Mpl 

<0RC(K)*5.0 
KX*  NPUL  S E 

IF(K.LT.NPULSE)KK=K 
DO  41=1, KK 
J*  K + l-l 

4 XORG(K)=XORGIK)+HPULSE< I ) *:VORG(  J ) 

5 CONTINUE 
30  TO  1303 

20  GO  TO  (101, 1"2,  103, 104, 105), IMRESP 

101  03  21K  = l . NPULSE 

21  HPJLSEtO-1.’* 

. GO  TO  2 

132  D3  22K  = 1 1 VPUL  SE 

22  HPJLSEU)=FLTAT(NPULSE*1-KI  /NPULSE 
GO  TO  2 

103  03  23K=1, NPULSE 

23  HPULSEI  K)=COS(  1.5 7.7796*FL04T( NPULSE M-H)  /NPULSE ) 

GO  TO  2 

104  03  24K=l, NPULSE 

24  HPJLSEIK)=EXPI-8.0*FL0AT( K*K-2*K+l)/ (NPULSE *NPULSE) ) 

GO  TO  2 

135  DO  25K=1, NPULSE 

25  HPULSE ( < ) = EXP ( -4.  0*FLOATI K-l ) /NPULSE ) 

GO  TO  2 

1030  RETURN 

END 
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SUBROUTINE: 
EQUATION : 


CORUPT 
AVGW  = 


AVGQ  = 


SUMQ2  = 


SUMW2  = 


SUMX2  = 


SUMV2  = 


_1 

XMP1 

_1 

XMP1 

J. 

XMP1 

_1 

XMP1 

_1 

XMP1 

_1 

XMP1 


MP1 

E 

k=l 

MP1 

E 

k=l 

MP1 

E 

k=l 

MP1 

E 

k=l 

MP1 

E 

k=l 

MP1 

E 

k=l 


WK 


QK 


QK 


WK 


x(k)' 


v(k)‘ 


Sample  mean  of  input  noise  seq. 


Sample  mean  of  output  noise  seq. 


Sample  mean  of  power  of  output  noise 


Sample  mean  power  of  input  noise 


Sample  mean  power  of  output  seq. 


Sample  mean  power  of  input  seq. 


SUBROUTINE:  CORUPT 

DESCRIPTION:  This  subroutine  calculates  the  standard  deviation  of  Input 
(SIGW)  noise  sequence  and  Output (SIGQ)  noise  sequence.  The 
values  of  the  mean  power  of  input  and  output  noise  along  with 
the  mean  power  of  input  and  output  sequence  are  also  calculated 


in  this  subroutine.  The  final  calculation  is  of  the  noise  to 


signal  power  ratio  of  both  input  and  output. 


PROGRAM  VARIABLES : MP1 

SIGQ 

SIGW 


V 

X 


NUMBER  OF  DATA  POINTS 

STANDARD  DEVIATION  OF  THF  OUTPUT  NOISE 
SEQUENCE  Q(k) 

STANDARD  DEVIATION  OF  THE  INPUT  NOISE 
SEQUENCE  W(k) 

CORRUPTED  INPUT  SEQUENCE 

CORRUPTED  OUTPUT  SEQUENCE 
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SUBROUTINE  CORUPT l X , V , S I GC, S I GW , HP  1 ) 

DIMENSION  <11 )> VII) 

RE AL*8  SUMV2, SUMX2,SUMQ2, SUMW2 , A VGQ , AVGW 

IX  INITIALIZES  UNIFORM  RAN00M  NUMBER  GENERATOR  IT8M  SUBROUTINE  RAMOUI 
RANDU  IS  GALLED  BY  SUBROUTINE  NORM  WHICH  GENERATES  NORMAL  DEVIATES 


INITIATE  T HE  RANDOM  SEQUENCE  GENERATOR 
IX=65549 
XMP 1 =MP I 
SUMX2  = 0.'>D0r? 

sjmv:=o.ddoo 

D039k=  i. » mpi 
SJMX2  = SJMX2*X (X  )*X(X ) 

SJMV2  = S JMV2+VIK  >*V(X> 

CONTINUE 

SUMV2=SUMV2/XMP1 
SUMX2= 3UMX2/XMP1 
IF(  SUMV2.EQ.''.GD0)SUMV3  = 1.0DD 
IF(  SUMX2.EQ.r'.  JD0)SUMX2  = l.(jD0 

FOR  INPUT  = 3,  S I GW  AND  SIGQ  BECOME  STD.  DEV.  CF  NOISE 

SI5W=DSQRT<  SUMV2*SIGW) 

SI3Q=DSQRTISUM<2*SIGQ) 

WRITE(6,897)SIGQ,SIGW 
897  FORMAT ( //30X, • SIGQ=" ,G17. 10 ,5X, • SIGH=« ,317. ID,//) 

SUMQ2=U.ODOC 
SUMW2»..0D30 
AVGQ=0. ODOO 
AVGW=0.0DCJ 
D340X=1,MP1 
CALL  NORMIWX.IX) 

CALL  NJRMIQK, IX) 

WK  = WK*S I GW 
QK=  QK*S  I GQ 
SJMQ2sSUM02+QK*QK 
SUMW2=SUMW2+WK*WK 
AV3W=A V3W+WK 
AVGQ=AVGQ+QK 
VI X ) =V ( X ) + WK 
40  X I X ) =X ( X ) + QK 

AV3W=AVGW/XMP1 

AV3Q=AVGQ/XM°l* 

SJMW2=SJMW2/YMP1 

SJMQ2=SJMQ2/YMP1 

ERR  X*DSQRT ( SUMQ2/SUMX2 ( *1GO.O 

ERRV=OSQRT(  SUM W 2/ SUM V 2 )*1*‘ O.G 

WRITE! &, 10U1) AVGQ.AVGW,SU''Q2,SUMW2.SUMX2,SUMV2.ERRX,ERRV 
1031  FORMATI///,2'*X, "SAMPLE  MEAN  OF  OUTPUT  NOISE  SEQUENCE  = a. Ell. At/ 


1,20X,"S 

AMPLE 

MEAN 

OF 

input 

NOT SE  SEQUENCE  = • , Ell. 4,/, 23X, 

2 " S AMPL  E 

MEAN 

POwER 

OF 

output 

NOISE 

",Ell.4,/,2JX, 

3 " SAMPLE 

MEAN 

POWER 

OF 

INPUT 

NOISE 

*tEll.4»/,2DXt 

A'SAMPLE 

MEAN 

POWER 

OF 

OUTPUT 

SEQUENCE  = 

• t £ 11 . 4 . / t 20 X, 

5" SAMPLE 

MEAN 

POWER 

OF 

input 

SEQUENCE  = 

" » E 11 . 4 1 / 1 2D  X, 

6*130.0  TIMES  THE  SQUARE  ROOT  OF  THE  NOISE  TO  SIGNAL  POWER  RATIO  OF 
7 THE  OUTPUT  = • ,F 7. 3 , / , 2C X. 

8 ’ 1 03 . 3 TIMES  THE  SQUARE  ROOT  OF  THE  NOISE  TO  SIGNAL  POWER  RATIO  OF 
9 THE  I NPUT  « • ,F7.3) 

RETURN 

END 
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SUBROUTINE : ERROR 

PURPOSE:  CALCULATE  PERCENT  MEAN  POWER  ERROR  IN  RECONSTRUCTION  AND 

PERCENT  OF  SQUARE  ROOT  OF  POWER  ERROR  IN  RECONSTRUCTION. 


EQUATIONS : 


AVGW  - I . (XORG(I)  - XREC(IH  2 * 
AVOW  L l X0RG(I)  J 

AVGQ  - VAVGW  * 100 
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SUBROUTINE:  ERROR 

DESCRIPTION:  The  subroutine  ERROR  calculates  the  Output  sequence  from  the 
Input  sequence  and  the  GAMMA  Matrix.  The  sequence  called 
XREC  is  compared  to  XORG,  the  original  output  sequence.  The 
sequence  XREC  is  generated  in  the  subroutine  RESPON.  The 
comparison  of  XREC  to  XORG  consists  of  calculating  the 
percent  mean  power  error  and  the  percent  of  square  root  of 
power  error. 


PROGRAM  VARIABLES : FDBACK 

GAMMA 

IDLY 

MP1 

N 

V 

XLAMDA 


VARIABLE  TO  PROVIDE  FEEDBACK  IF  DESIRED 

MEASUREMENT  VECTOR 

DELAY  INTRODUCED  IN  INPUT  NUMERATOR 

NUMBER  OF  DATA  POINTS 

ORDER  OF  SYSTEM 

CORRUPTED  INPUT  SEQUENCE 

WORKING  ARRAY 


XREC 


RECONSTRUCTED  OUTPUT  SEQUENCE 
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SUBROUTINE  ERROR! XREC.V, GAMMA ,MP1 ,N, XLAMDA, XORG .FDBACK, IDLY) 
DIMENSION  XR'Cll) ,V( 1 ) ,XORG(l  ) 

DIMENSION  VVVI20) 

REAL*B  GAMMA! 1)  , XL AMD All) ,AVGW, SUMV2  » AVGO 
INTEGER  FDBACK 

CALL  RESPON(XREC»V,N, GAMMA, XL AMDA, MP 1 , FDBAC K I 
AVGW*0. 3000 
SUMV2*3#  ODD 
00261=1, MP1 

S'JMV2*SJMV2+<ORG(  I)*XORG(  II 
AVGQ*XORG< I )-XREC( I I 
4VGW*AVGW+AVGQ*AVGQ 
AVGW=AVGW/SUMV2 
AVGC»DSORT (AVGW  I 
AVGQ=103.0*AVGQ 
AV3W=1C0.3*AVGW 
WRITEI6,27)AVGW,AVGQ 

FORMAT (IX,*  PCR  CENT  MEAN  POWER  ERROR  OF  RECONSTRUCT  I ON* , F8. 3, /// , 
11X, • PER  CENT  OF  SQUARE  ROOT  OF  POWER  ERROR  IN  RECOSTRUCT ION* , F8. 3) 
RETURN 
ENO 
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SUBROUTINE: 

PURPOSE: 


FLOW  CHART: 


FILLV 

GENERATION  OF  DISCRETE  POINTS  FOR  A VARIETY  OF  WAVEFORMS 
(For  the  correspondence  of  the  waveshapes  and  input  parameter 
INPT  see  page  ) 
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SUBROUTINE : FILLV 


8)  TRIANGULAR  WAVE 
PERIOD  = NPUL 


9)  DECAYING  EXP  + 

DAMPED  SIMUSOID 

T = NPUL 
1 

T = 2*NPUL, 

2 

PERIOD=l . 91*NPUL 


11),  12),  13)lNDEPEN)ENT| 
RANDOM  NOISE 
GENERATORS 


IF  ISHIFT  = 1 
WAVEFORM  SHIFTS  NPUL  T0| 
THE  RIGHT  AND  ZERO 
REPLACES  FIRST  NPUL 

IQIH1S 


DESCRIPTION:  This  subroutine  builds  an  array  of  NPT  points  defined  by  the 
choosen  wave  form  and  parameter  (NPUL)  of  that  waveform.  It 
is  useful  in  approximating  input  signals  for  excitation  of 


control 

system. 

PROGRAM  VARIABLES: 

INPUT 

DESIRED  WAVEFORM  OPTION 

NPT 

NUMBER  OF  DATA  POINTS 

NPUL 

WAVEFORM  PARAMETER 

V 

GENERATED  OUTPUT  SEQUENCE 
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SUBROUTINE: 


FILLV 


1)  UNIT  PULSE 

2)  UNIT  STEP 

1 

3)  DOUBLET 

-*l 

WIDTH  = NPUL  , 

1 

\h 

4)  SQUARE  WAVE 

— *1  NPUL 

PERIOD=NPUL  1 

iruir 

5)  SQUARE  WAVE  TO 

EXPONENTIAL  DECAY  1 

— NPUL 

1 n n- 

PERIOD  = NPUL 

TIME  CONSTANT=NPUL 

u U 

6)  EXPONENTIAL 

TIME  CONSTANT =NPUL 

k. 

7)  PERIODIC  IMPULSE 
PERIOD=NPUL 

1— *1  NPUL 

t T t T 

8)  TRIANGULAR  WAVE 
PERIOD=NPUL  , 

I* *1 

A A j 

1 V V 

9)  DECAYING  EXP  + 
DAMPED  SINUSOID 

T i = NPUL 

YlpPRflklMATE 

= 2*NPUL,  Period 

= I . 91*NPUL) 
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SUBROUTINE  FILL  V(  V,  NP  T , I NPUT,  NPUL ) 

C FILLS  THE  ARRAV  FOR  INPUT  ACCORDING  TO  INPUT  OPTION  DESIGNATED 

DIMENSION  VU  ) 

GO  TO  (1.2.3, A, 5, 6, 7,8,9, iU), INPUT 

1 mi>u 

03  101  I *2 ,N°T 
191  (1110.3 

GO  TO  999 

2 03  102  I * 1 , NP T 

102  V( I 1*1.3 
GO  TO  999 

3 DO  103  I-l.NPT 
V( I 1-3.0 

IF( I .LT.NPUL /2 ) V( I ) -1 • 

IF ( I .GE.NPUL. AND. I «LT • NPUL)  V(I)--1.0 

103  CONTINUE 
GO  TO  999 

A 03  104  I-l.NPT 

TIII'1.9 
TPUL-I /NPUL 

I F ( I-TPUL*NPUL. GE .NPUL/2 I VIII— 1.0 

104  CONTINUE 
GO  TO  999 

5 00  105  l-l.NPT 
V( I ) = 1.9 

1F(  I.Gc. NPUL/2. ANO.I. LT.NPUL)  Vdt  — 1.0 
I FI  I.GE.  1 1.5)  *NPUL.  ANO.I.  LT.2*NPUU  V(I)— 1.0 
ARG1-2.5-FL0&TI I ) /FLOAT ( NPUL ) 

IF< I.Gc. <2.51 *NPUL)  V< I)=EXP< ARG1) 

105  CONTINUE 
GO  TO  999 

6 03  136  I-l.NPT 

ARG2  — FLOAT  1 1 ) /FLOAT  I NPUL  I 

106  VII >=EXP(ARG?) 

GO  TO  999 

7 DO  107  I *1 ,N°T 

107  vm-o.o 
N-NPT/NPUL 

vrn-i.3 

03  1G71  J=1,N 
I-J*NPUL 
1971  VI I )-1.9 

GO  TO  999 

8 03  108  I-l.NPT 
TPJL-I/NPUL 
ITPUL-TPUL 

VU )«( 2. *FLOAT( D/FLOAT  I NPUL )-2.*TPUL)*( -1.0) **ITPUL 
IF! I-TPUL*NPUL.GE. NPUL/2)  V(I )«2*(1*TPUL-FL04T( 1 ) /FLOAT ( NPUL) >*( 
1.9 ) **I TPUL 

108  CONTINUE 
G3  TO  999 

9 00  109  I-l.NPT 

ARG3— FLOAT  ( I ) / FLOAT  ( NPUL  ) 

ARG4»“.5*FL0AT( I) /FLOAT (NPUL) 

ARG5-3.296*FL0AT(I ) /FLOAT  INPUL ) 


IMPULSE 

STEP 

DOUBLET 

SO  WAVE 

SQ-EXP 

EXP 

PRO  IMPL 

TRI  HAVE 

1 

EXP+OSC 
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13V  V ( l )=E<P(ARG*I«S*P<ARG4)*SIN(ARG5> 
G3  TO  999 

13  (<*629327213 

03  113  1=1, NPT 

4 = 3. 

03  1101  K=l,’2 

I Y= I <*65539 

I F ( IY)  1102, 1133,11  .'3 

1132  IY=  IY«2U7493667*1 

1133  YF L = I Y 

YFL=YFL«. 465661 3£-9 

1<=IY 

A=A*YFL 

1131  CONTINUE 

Yi I ) = A-S.O 
113  CONTINUE 
53  TO  999 
999  CONTINUE 
RETURN 
END 
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SUBROUTINE:  FINDQ 

PURPOSE:  CALCULATES  MEASUREMENT  FILTER  POLE(Q)  AND  NUMERATOR  (DEL) 

EQUATION:  H (Z)  = — ^7 

m 1-QZ'1 

FLOW  CHART: 


NCSL  TM-204-77 


SUBROUTINE:  FINDO 

DESCRIPTION:  FINDQ  determines  the  measurement  filter  pole  and  numerator. 

The  subroutine  uses  an  iterative  process  of  calculating  DEL 
and  Q.  The  iteration  is  satisfied  when  the  variable  POW  is 
within  + .5%  TEM. 

PROGRAM  VARIABLES:  DEL  NUMERATOR  OF  1st  ORDER  MEASUREMENT 

DIGITAL  FILTER 

MPl  NUMBER  OF  DATA  POINTS 

N ORDER  OF  MODEL 

Q 


XI 


MEASUREMENT  FILTER  POLE 
COEFFICIENT  VECTOR  (same  as  Gamma) 
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SUBROUTINE  FtNDQI Q.OEL , XI , X , MP1 , N) 

OMENS  ION  X<1) 

REAL*B  9(11  ,PELm.TEH,POH,PT,QS,OBIG,aSKAL,Xl(l),SUM 
SUH-0. 300 
D0667K-1 » NP1 
667  SUN- SUN»X(K>**2 

NPl -N* l 
0311-1. N 
LP1-L+1 

TEN-iO'J.ODO/DFLOAT  l NP1 ) •DFLOATI NP1-L) 

BBIG-1.0D0 

qsmal-;-.30g 

100  0S»(QBl5*0SMAU/2.rD0 
QILI-OS 

DEL ( L ) -1 . 000-QS 
PT -0.000 
D04I-1.LP1- 
6 XKII-O.  ?DO 

D03K-1.MP1 
Xll  ll-XUJ 
D35I-1.L 

5 XII  I*1)-X1II+1)*QII)+X1II  )*DEUI> 

3 PT»PT4Xl(LPn*Xl(LPl) 

P3N-PT/SUN*lr,0. 000 

IF(P0W.LE.l.n35D3*reN.AN£>.POW.GE.7eN*O.99500)SO  TO  1 
I F ( POM. GT .TEN) GO  TO  6 
OBIG-OS 
GO  TO  10!) 

6 QSNAl-QS 
GO  TO  103 

1 CONTINUE 

RETURN 
ENO 
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SUBROUTINE:  G RAMI I 

PURPOSE:  PERFORMS  GRAM  II  TECHNIQUE  WHICH  YIELDS  THE  GRAM  MATRIX  (G) , 

NOISE  CORRECTION  MATRIX  (Z)  AND  THE  COEFFICIENT  MATRIX  (GAMMA) . 

GAMMAt (k ) *GAMMAj (k ) 


EQUATION: 


MP1 

G,  . = Z 
ij  k=! 


-A26- 


A 


r 


I 


NCSL  I'M- 204- 7 7 

SUBROUTINE  : GRAMII 
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SUBROUTINE:  GRAMII 
PROGRAM  VARIABLES: 

DELTA 

GAMMA 

G 

IDLY 

ILEVIN 

IZTS 

MAX 

N 

QOPT 

QSAV 

RHO 

SIQQ 

SIWW 


SAMPLING  INTERVAL 
MEASUREMENT  VECTOR 
G MATRIX 

DELAY  INTRODUCTED  ON  INPUT  NUMERATOR 

VALUE  IS  EITHER  0 OR  1.  0 GRAM 

TECHNIQUE  IS  PERFORMED.  1 LEVIN 
TECHNIQUE  PERFORMED. 

SEE  SUBROUTINE  ZTOS 

DIMENSION  SIZE  OF  2 DIM  ARRAYS  IN 
THE  DIMENSION  STATEMENT. 

ORDER  OF  SYSTEM 

MEASUREMENT  FILTER  OPTION.  QOPT=l  Q(I) 
GENERATED  IN  FINDQ  QOPT=0  Q(I)=QSAV 

MEASUREMENT  FILTER  POLE 

EXPECTATION  OF  (W(K)*Q(K)) 

STANDARD  DEVIATION  OF  OUTPUT  NOISE  SEQUENCE 

STANDARD  DEVIATION  OF  INPUT  NOISE  SEQUENCE 
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SUBROUTINE:  GRAMII 

PROGRAM  VARIABLES  CONTINUED: 


V CORRUPTED  INPUT  SEQUENCE 

X CORRUPTED  OUTPUT  SEQUENCE 

XLAMDA  EIGENVECTOR 

NOISE  CORRECTION  MATRIX 


Z 


1 
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SUBROUTINE:  GRAMII 


I 


CALCULATION  OF  GRAM  MATRIX  "G" 
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SUBROUTINE  GRAM  1 1 1 X, V, MP l , S ICQ, S IWW.RHO.N, DELTA .QSAV.QOPT , I REM, 
l JZTS, GAMMA, XL AMOA.G.Z, MAX, ILEVIN, IDLY, INORM) 

C 

C THIS  SUBROUTINE  PERFORMS  THE  GRAM  II  TECHNIQUE 
C 

DIMENSION  XU)  , V<1  >,  G(  MA  X.l),  Z!  MAX,  1J,  GAMMA  II),  XLAMDAU  >,0(201, 
10EL  < 20 1 

DOUBLE  PRECISION  G, Z .GAMMA, XLAMDA.OELTA, OEL ,PR00,0 ,QSAV 
INTEGER  QOPT 

REAL*8  XMEAN,  GAM! 25 ) « FAC 
REAL*8  S,F,S1,S2,VARQ, VARW.GI 

COMMON  /MATRIX/S!20.20>,F!20,20>,GI(20,20J.S1!10,10> ,S2(10,10) 

MAX2=MAX/2 

WRITE! 6 « 1000) 

1000  FORMAT! 1H1,20X.'THE  GRAM  II  TECHNIQUE') 

C JOPT  » 0.  IF  DIRECT  TRANSMISSION  IS  ASSUMMED 

J0PT«0 

C JOPT  * l IF  NO  OIRECT  TRANSMISSION  IS  ASSUMMEO 

IF! IREM.NE.O>JOPT»l 

C 

C IOPT  » 0 NOISE  ON  BOTH  INPUT  ANO  OUTPUT  IS  ASSUMEO 

C IOPT  * 1 NOISE  ON  OUTPUT  ONLV  IS  ASSUMED 

C IOPT  - 2 NOISE  ON  INPUT  ONLV  IS  ASSUMEO 

C 

IOPT*0 

IFISIWW.EQ.O.OHOPT-l 

IFISIQQ.EQ. 0.0. AND. SIMM. NE. 0.0 ) I CPT«2 

C 

C OEL  IS  THE  NUMERATOR  OF  THE  KNOWN  FIRST  ORDER  DIGITAL  FILTERS 

IFIQOPT.NE.OI  GO  TO  21 
00191-1. N 

DEL  I I ) «1. OOON-QSAV 
19  Qi I I-QSAV 
GO  TO  22 

21  CALL  FINOQIQ, DEL, GAMMA. X.MPl.N) 

22  CONTINUE 

WRITE!*. 2020)  . 

2020  FORMAT! lOX.'Q  PARAMETERS' ) 

CALL  PRVECIQ.N) 

NPl -N-l 
NP2«N»2 
NPNP2«N*N^2 
NR-NP1-IR6M 
NP1PIR-NP1MREM 
VAR W> 0.0 
VARQ-0.0 
003001-1 »MP1 
VARW-VARW*V! I l*V! 1 1 
300  VAAQ-VARQ»XII)*XI1> 

VAR Q*0 SORT  I VAR0/MP1 1 
VARW-0SQRTIVARW/MP1I 
SI GQ-S1 QO/VARQ 
SIGW«SIWW/VARW 

IF! SIWW.EQ.O.O. ANO.SIQQ.EQ.O.OI SIGO-l.O 
NPNP2-NPNP2*! 


-A31- 


NCSL  TM-204-77 


NR=NR+l 
00101*1. NPNP? 

GAM!  11=0.000 
GAMMA! I 1*0.0000 
0010 J=  I , NPNP? 

10  G( I . J 1 =0.0000 
GAHl 11*1.000 
C 

C CALCULATING  THE  G HATR IX 

D050K* 1 , MP 1 

IF! ILEVIN.EQ. 11G0  TO  31 
IF!K-IDLY)25,25,26 

25  GAHMA!NP2 1*0.0000 
GO  TO  26 

26  F AC* 1. 0 

IF! INORM.EQ.l )FAC=VARW 
GAMMA !NP21*V!K- IDLY) /FAC 
FAC=1.0 

IF! INORM.EQ.l )FAC*VARQ 
GAMMA! 1 ]*X!K)/FAC 
26  CONTINUE 

00301=1, N 

GAM  ( I ♦ 1 1 *GAM(  1*1 1*0!  I 1 ♦GAM!  I )*OEHI) 

GAMMA! I*l)*GAMMA!I >*DELt I )♦ GAMMA  I I ♦ 1 1 *Q!  I ) 

30  GAMMA! I*NP2)  = GAMMA!I+NP1)*0EL!I  1 ♦GAMMA! I+NP2)*0( 1 1 
GO  TO  35 

31  CONTINUE 

00  32  1*1, NP1 
GAM! 1+1)=GAM! I) 

IFIK-IDLY-I.LT.O)  GAMMA! I 1*0.0000 
IF! K-IOLY-I.LT.O)  GAMMA! I +NP1 1*0.0000 
IFIK-IDLY-I.LT.O)  GO  TO  32 
FAC=1.0 

IF! INORM.EQ.l )FAC=VARQ 
GAMMA! 1)=X!K* 1-IJ /FAC 
FAC=1.0 

IF! INORM.EQ.l )FAC*VARW 
GAMMA! I+NP1)*V(K+1-IDLY-I  l/FAC 

32  CONTINUE 

35  CONTINUE 

GAMMA! NPNP2)*GAM!NP1) 

D060  I = 1»NPNP2 
0060 J= I , NPNP 2 

60  G < I • J I *GI 1 , J) ♦GAMMA! I I 'GAMMA! J) 

50  CONTINUE 

00601=2, NPNP? 

K*I-1 
0060 J* 1 , K 

60  G!  I ,JJ  = G!  J,  I) 

HRITE16, 10021 

1002  F0RMAT(20X,'THE  G MATRIX*  1 

CALL  PRMAT!G,NPNP2,NPNP2,MAX) 

C 

C CALCULATING  THE  NOISE  CORRECTION  MATRIX  l BY  SUBROUTINE  BUILOZ 

CALL  BUIL02!Z,S,GAMMA,N,MP1,SIGW,SIGC,RH0,0EL,Q,MAX,ILEVIN1 
IF!  JOPT)70,<)0,70 
70  CONTINUE 

IF! ILEV IN. EQ. 1 1 GO  TO  81 
D0B0J*1,NPNP? 

00801*1, NR 

ZINPlM,J)*Z!NPiPIR»I«J) 
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60  G!NP1*I,JI-G!NPIPIR*I,J) 

NPNP2-NPNP2-IREM 
D035J-1.NPNP2 
00851*1, NR 

ZIJ,NP1*I)»ZI  J.NP1PIRM) 

85  G1 J.NP1*I)*G! J,NP1PIR*I) 

GO  TO  90 

61  00  82  J-1.NPNP2 
ZINP1*NR*1, J)«G!NPNP2, J1 

82  G!NP1*NR*1,  J)*GINPNP2,  J> 

NPNP2*NPNP2-!REH 
00831*1, NPNP2 

Z! I » NP 1 *NR* 1 1 * Z 1 1 , NP 1*NP 1 *1 ) 

83  G(I,NP1*NR*1)-GII,NP1*NP1*1> 

CALL  PRHAT < G, NPNP2 ,NPNP2 , PAX  I 

90  NPNP2-NPNP2-1 

IF  < 10PT-1 1617 ,605,618 
C 

C NOISE  ON  BOTH  INPUT  ANO  OUTPUT 

617  CALL  S0LVE21Z.G, GAMMA, XLAPDA.NPNP2, 1 . MAX) 
XME AN-XL AMDAI NPNP2+1 ) 

NR*NR-1 
GO  TO  6C6 
C 

C NOISE  ON  OUTPUT  ONLY 

605  CALL  S0LVE21Z.G, GAMMA, XLAM0A,NP1, NR, PAX) 
XMEAN-XLAMDAINPNP2*1) 

NR-NR-l 
GO  TO  606 
C 

C NOISE  ON  INPUT  ONLY 

618  NPP»NPNP2*i 
NR*NR-1 

005 501*1, NPP 

550  GAMMA! I l-GI I, NPP) 

D0551J-1.NR 

JJ-NR-J*! 

D0551l-1,NPP 

551  Gll,NP2*JJ)*GII,NPl*JJ) 

D0552I*l,NPP 

552  G ( I , NP2 ) -GAMMA! I ) 

005531*1, NPP 

553  GAMMA(I)*GINPP,I) 

005 541-1, NR 

II-NR-I*1 

D0554J-1.NPP 

554  GINP2+II.J)-G!NPi*II.J) 

005551-1, NPP 

555  GI NP2, I )- GAMMA! I ) 

CALL  S0LVE3! Z,G, GAMMA, XLAM0A,NP2, NR, PAX) 
XMEAN-XLAMDAINP2) 

D0556I -1,NR 

556  XLAMDAI NP1* 1) -XLAMOA! NP2* I ) 

606  IF! JOPT  >120,130,120 
120  NPNP2-NPNP2MREM 

IF!  ILEVIN.EQ.DGO  TO  124 
001221-1, NR 

122  XLAMOA ( NPNP2-I*1 l-XLAMOA! NP2*NR-I ) 
001231-1, IREM 

123  XLAMOA! NP 1*1 ) -0,0000 
GO  TO  130 


r 
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124  NN3MIR=NPNP2*l-IREM 
00I25l=NN3MIR,NPNP2 

125  XLAMOA  t I)=0«n000 
130  CONTINUE 

FAC-1.0 

IF!  INORM.EQ.l )FAC=VARQ/VARW 
00301t'NP2,NPNP2 
301  XL A MOA ( I )=XLAMDA! I I *FAC 
WRITEI6.1001) 

1001  F0RMAT!10X,*TME  synthetic  coefficient  vector,  xlamoa,  ism 
CALL  PRVEC1XLAMOA.NPNP2) 

DO  15C  I=1,NPNP2 
150  GAMMA! I )* XLAMOA! 1 1 

IF( ILEVIN.EQ.11C0  TO  165 
C 

C GENERATING  GAMMA  FROM  XLAMOA 

CALL  BUIL0A(S,Q,0EL,N,MAX» 

D0160I  = 1 « NPNp2 
GAMMAC 11=0.0000 
FAC*1.0 

IFF INORM. EO. 1 )FAC=VARO 
00160J=1,NPNP2 

160  GAMMA! I ) = GAMMA( I) ♦ S ( 1 . J ) * XL AMOA < J ( 

XMEAN=XMEAN*S! 1 ,NP 1 > *FAC /GAMMA! 1 ) 

165  CONTINUE 

WRITE16,655)XMEAN 

655  FORMAT!/lX,*MEAN  COEFFICIENT  IS  *,013.6,//) 

D0200I *2, NPNP2 

200  GAMMA! n=GAMMA( I)/GAMHAI1 I 
GAMMA! 1 1*1.0000 
IF! ICLY.EQ.OtGO  TO  172 
IDLY1* I0LY*1 

00  170  II=IDLVl,N°l 

1 = NPfiP2*l-I  I 

170  GAMMA! I ♦ I OL Y) =G AMM A! I ) 

00  172  1*1, IOL Y 
GAMMA! I+NP1) =0.0000 
172  CONTINUE 
C 

C CALCULATING  THE  EQUIVALENT  CONTINUOUS  DESCRIPTION 

CALL  I ZTOS! GAMMA.N , OELTA * IZTS ) 

WRITE! 6,1003) 

1003  FORMAT!///, IX, 1QP11H-) , / , IX , ICO  1 1H-) ) 

RETURN 

END 
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SUBROUTINE: 

PURPOSE: 

EQUATION: 

FLOW  CHART: 


IZTOS 

SEPARATES  THE  NUMERATOR  FROM  THE  DENOMINATOR  PARAMETERS  IN 
GAMMA. 

DENOMINATOR;  XI (I)  = GAMMA(I) 

NUMERATOR;  X2(I)  = GAMMA (NP1  + I) 
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SUBROUTINE:  IZTOS 

DESCRIPTION:  This  subroutine  takes  the  coefficient  vector  GAMMA  and 

separates  the  vector  into  the  numerator  (X2(I)  array)  and 
denominator  (X1(I)  array). 

PROGRAM  VARIANCE:  DELTA  SAMPLING  INTERVAL 

GAMMA  COEFFICIENT  VECTOR 

IZTS  « 0 Z DOMAIN  TO  S DOMAIN  CONVERSION  NOT 

PERFORMED 

* 1 LOGARITHMIC  TRANSFORMATION  IS  PERFORMED 

= 2 PULSE  DELAYED  TRANSFORMATION  IS  PERFORMED 

N SYSTEM  ORDER 
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SUBROUTINE  I2TOS(GAMMA,N, DELTA,  IZTSI 

I2TOS  SEPARATES  THE  NUMERATOR  FROM  THE  DENOMINATOR  PARAMETERS 
IN  GAMMA 

DIMENSION  GAMHA(l»,xmO)  ,X2(10) 

DOUBLE  PRECISION  GAMMA, XI ,X2, DELTA 
NP1»N*1 
200  D03 I -1 , NP1 

Xl<  n = GAMMAU) 

3 X2(I1»-GAMMAINP1*I » 

CALL  2T0SCX1,X2,N,0ELTA, IZTS) 

IZTS*I2TS*1 

IF! 12TS.E0.4t  GO  TO  200 

RETURN 

END 
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SUBROUTINE:  POLCON 

PURPOSE:  CONSTRUCTS  POLYNOMIAL  FROM  ITS  ROOTS 

FLOW  CHART: 


DESCRIPTION:  This  subroutine  constructs  a polynomial  from  its  roots  and 

the  polynomial  coefficients  are  stored  in  an  array  R2(I) 
in  ascending  order. 

PROGRAM  VARIABLES  C ROOTS  USED  TO  FORM  POLYNOMIAL 

K OPTION  WHICH  SUPPRESSES  POLYNOMIAL  CONSTRUCTION 

WITH  SPECIFIED  ROOT(S) 

N ORDER  OF  SYSTEM 


R2 


COEFFICIENTS  OF  POLYNOMIAL 
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SUBROUTINE  P0LCGN(C,R2,K,N( 

A POLYNOMIAL  CONSTRUCT 1 ON  PROGRAM  NEEDED  FOR  ZTOS 

DIMENSION  CI1),R2(1) 

C3MPL  E A* 1 6 C.R2.COMP 
REAL«8  DEI?) 

EQUIVALENCE  (COMP, DC) 

NPl=N+i 
D31 01=2 , NP1 
R2(  I ) = '..0D r: 

R2  1 1 ) = 1 . 5D:>.‘ 

D3A I = 1 , N 
C3MP  = C(  I ) 

IF(  I.EO.K.OR.  (DC(I).EO.'J.'  DC1.  ANO.  DC  l 2 ) . EO.  3.  JDD)  )GO  TO  A 

D02 J J=1 , I 

J=I-JJ*1 

R2(J*l)=R2U+l)«Cm+Rc<  J> 

R2 ( 1 ) = R2 ( 1 )*CI I ) 

CONTINUE 

RETURN 

END 
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SUBROUTINE:  PRCVEC 

PURPOSE:  This  subroutine  prints  out  a complex  single  dimensioned  array. 

EQUATION:  Complex  number  A + BJ  is  printed  (A,  BJ) 

FLOW  CHART: 


DESCRIPTION:  This  subroutine  is  called  in  the  ZTOS  subroutine  when  the 

poles  and  zeroes  in  the  S domain  are  needed. 

PROGRAM  VARIABLES:  A ARRAY  TO  BE  OUTPUT 

N NUMBER  OF  ELEMENT  IN  ARRAY 


NCSL  TM- 204-77 


- - — m ^ 


I 


SUBROUTINE  P«CVECIS,N> 

*• 

l THIS  SUBROUTINE  PRINTS  OUT  A COMPLEX  SINGLE  DIMENSIONED  ARRAY 

Z A COMPLEX  NUMBER  OF  THE  FORM  A ♦ B J IS  OUTPUTTED  IN  THE  FORM 

; I A,  M J)  WHcRc  J = SQUARE  ROOT  OF  -1 

DIMENSION  All  ) 

COMPLEX*!*  A 
WRITE!*, 2) 

WRITE46.1) 1 AI  I)  * I = 1 * N ) 

1 F0RMAT(1X,1H( ,D17. 10. 1H, ,D17.1D,3H  J ) ) 

WRITEI6.2) 

MR  I TCI  * , 2 ) 

2 FORMAT  1 1 ) 

RETURN 

END 
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SUBROUTINE:  PFMAT 

PURPOSE:  SUBROUTINE  OUTPUTS  DOUBLE  PRECISION  DOUBLE  DIMENSION  ARRAY 

FLOW  CHART: 


DESCRIPTION:  This  subroutine  is  called  in  GRAM  II  and  takes  an 

array  of  two  dimensions  and  gives  an  output  of  the  same  two 
dimensional  array  in  double  precision. 


PROGRAM  VARIABLES:  A 

M 

N 

NMAX 


OUTPUT  DOUBLE  PRECISION  ARRAY 

MATRIX  "A"  COLUMN  DIMENSION 

MATRIX  "A"  ROW  DIMENSION 

DIMENSION  SIZE  OF  TWO  DIMENSIONAL  ARRAY 
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SUBROUTINE  PRMAT(A,N,M,NMA<) 

DOUBLE  PRECISION  A 

THIS  SUBROUTINE  OUTPUTS  DOUBLE  PRECISION  DOUBLE  DIMENSIONED  ARRAY 
DIMENSION  A INMAX* II 
URlTE(6,l> 

002  I *1 1 N 

2 WAITEIo.3) <A( I, 

3 FORMAT ( IX. 10013.5 » 

HR  I TE ( 6> 1 ) 

HRITE(fa.l) 

1 FORMAT!/) 

RETURN 

END 
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SUBROUTINE:  PRVEC 

PURPOSE:  SUBROUTINE  OUTPUTS  DOUBLE  PRECISION  SINGLE  DIMENSION  ARRAY 

FLOW  CHART: 


DESCRIPTION:  This  subroutine  is  called  in  GRAMn  and  other  routines  to  print 
one-dimensional  arrays  in  double  precision. 

PROGRAM  VARIABLES:  A ARRAY  THAT  IS  OUTPUT  IN  DOUBLE  PRECISION 

N NUMBER  OF  ELEMENTS  IN  ARRAY 
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SUBROUTINE  P°VECtA.N) 

THIS  SUBROUTINE  OUTPUTS  DOUBLE  PRECISION  SINSLE  DIMENSIONED  AR«*V 

DIMENSION  At1  ) 

DOUBLE  PRECISION  a 
HRITttL.il> 

HRlTEtu.il I At I> i I* 1 *N ) 

FORMiTt K. 10013. 5) 

WRITE- 16.311 
WR  I T c I 6 « 31 1 
FORMAT  t / > 

RETURN 

END 


r 


SUBROUTINE: 

PURPOSE: 

EQUATIONS : 


NCSL  TM-204-77 

RESPON 

CALCULATES  RESPONSE  (X^  FROM  COEFFICIENT  VECTOR  (GAMMA) 
TIMES  THE  ARRAY  XLAMDA. 

[Xv,  (XLAMDA)]  [ GAMMA  ] = 0 

In. 


or 


XK  « - tVr-'Vn  Y • • Vn1 


Y(2) 


|_y(N+N+22| 


also  XK  = - (XLAMDA)  (GAMMA) 


FLOW  CHART: 


/V 

✓ V 


*<eOMPS(2)-0.j>-N- 

N ✓ 

^ / 

V 

i 

_+_Yes 

COMPs""(2)“I .0 1 

• COMPS(1)-0.0  I 

> | 

l C0MPS(3)-0.0  . 

' . J 


INITIALIZE  XLAMDA 
XSAV 


XLAMDA(J)  - XLAMDA(J-l) 
XLAMDA (1)  - XSAV 
XLAMDA(NPl)  “ V (K) 


1 

♦ 


* 

r w 1 

I INITIALIZE  4>  , 


XSAV  - 0.0  

XSAV-XSAV- (XLAMDA) T (GAMMA 


XSAV 

EKM1 

WKM1 


A — 


Yes 


*- 1 J 

I 

X 

X \ 

-<FDBACK,  LT . 1 ^ 

v : " 


T No 


[ek 


V (K)  - GAIN  * XSAV  1 
WK*COMP  S ( I ) *WKM1+C0MP  S ( 2 ) I 
-COMP  S ( 3 ) *EKM1  | 


| EKM1 
I WKM1 


(V(K)_ 


EK 

WK 

WK 
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SUBROUTINE: 

DESCRIPTION: 


RESPON 


Subroutine  RESPON  determines  the  response  of 

-n 


H(z)  - 


b + b.z  ^ + b_z  ^ + ....  + b z 
o 1 2 n 

. . “1  . ~2  . . -n 

1 + a, z + a,z  + ....  + a z 
12  n 


to  the 


input  sequence  V(K).  The  coeffieints  are  entered  as  an  NPNP2  ■ 


N + N + 2 vector  CAMMA 
PROGRAM  VARIABLES:  FDBACK 

GAMMA 

MP1 

N 

V 

X 

XLAMDA 


- ( 1 • a j aR,  - bQ "bp) 

NO  FEEDBACK  FDBACK  - 0 
NEGATIVE  FEEDBACK  FDBACK  - 1 

COEFFICIENT  VECTOR 

NUMBER  OF  DATA  POINTS 

MODEL  ORDER 

INPUT  SEQUENCE 

OUTPUT  SEQUENCE 

WORKING  ARRAY 


FEEDBACK  AND  COMPENSATION: 
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SUBROUTINE:  RESPON 

PURPOSE:  CALCULATES  RESPONSE  (X^)  FOR  THE  SYSTEM  SHOWN  IN  ABOVE  FIGURE. 

THIS  ADDITION  TO  SUBROUTINE  RESPON  INCORPORATES  THE  El  EXAB I L 1 TY 
OF  ADDING  NEGATIVE  FEEDBACK  AND  CASCADE  COMPENSATOR  IN  THE 
FORWARD  LOOP  FOR  OPTIHUATION. 

EQUATION:  EK  - - AJ^ 

WK  - WKMl  * CCJMFS  (I)  ♦ EK  * COMPS  (2)  - COWS  (S)*KRMl 

FLOW  CHART:  (Saa  flow  chart  for  RESPON  tha  dottad  •action  la  for  Eaadha.  » 
and  Cowpanaatlon  natworh.) 


r 


NCSL  W-204-77 


BEST  AVAILABLE  COPY 


Lj«  jllil  M|#OH  4,4,  I 

OI"«  .»  )*<  It  * l . 41 ; i .&*-«» I L t , lc 4*341 1 1 
M»l»*  l\44,*.4*«4,  IC.4*  4.G4IN 
hltill  »0*4'4 

t«|  I CMtl  »94  4 IVMil  MIMM  »»!»!• 

•110*4.4  j4  I S C Hi  *0  '641. • II  0tV|**MltD 


4|4l »4  t’4»» 

;j*«  > i 

,4IN**«II 

'i* : • 

44|  >M| 

lltlKO 


4100  HI  l,8WV0l  III,  800411131  .(0*  III  I 


<••1  *»*i 
l»t»V  . 

4*40 l • 

i.  4*  <4t  i • * 

t • • c *•»♦! n.M.  *• 

JHUf 

c »»•#.* 

il4fl,NN 

•••l*t. 

•••I*!. 

0J’ I 

|flMi:itiMIM  »0  111 

M«*|4  |.„4l4*lt44 

•■•c  • i i •**••.  ♦cD**im*i4-co**»i 

M4|  • 

f««l*|4 

414 

331  1 1 • 1 •'•<*1 
J*«»l*l 

4C  4 <104  I J ••  U4».UU*  l I 

03fll*i,4 

4c4*t>4l  J(»ul»34|  j.JI 
4144C4I 1 ••4$4V 
H4MU4I  4*1  1*414  1 
4S44*0.3033 
03*»l»l.4W»01 

<S4V*<S44-G4«*AI  | « 1 IML.4N04I  I I 
IM  04#SI  4|44).GT. 1,001  MOO  VO  IT 
<U  1*1144 
MTURN 


002  8 1 *4*  40| 
<111*0.0 
RETURN 
END 


' .4 
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SUBROUTINE: 

PURPOSE: 

FLOW  CHART: 


SOLVE1 

FINDS  MAXIMUM  EIGENVALUE  AND  CORRESPONDING  EIGENVECTOR  OF 
(Z-E*C)*V-0 


SUBROUTINE:  SOLVE1 
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DESCRIPTION:  Subroutine  SOLVE1  uses  the  first  quadrant  of  the  GRAM  matrix 

and  noise  correction  matrix  " " to  calculate  the  output  portion  of 
the  EIGENVECTOR  (N+l  elements)  plus  the  maximum  EIGENVALUE. 


PROGRAM  VARIABLES: 


G 

MAX 

N 

S21 

V 

Z 


FIRST  QUADRANT  OF  GRAM  MATRIX 
MAXIMUM  ROWS  PERMISSIBLE 
N+l  (ORDER  OF  SYSTEM  + 1) 
COEFFICIENT  MATRIX 
EIGENVECTOR 

NOISE  CORRECTION  MATRIX 
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SUBROUTINE  S*LVE1U,G,S21,V,N,MAXI 

FINOS  MAXIMUM  EIGENVALUE  AND  CORRESPONDING  EIGENVECTOR  OF 
I 2 - E*G  ) * V * C 
WHERE  l AND  G ARE  M X N MATRICES 

E IS  THE  e I GENVALUE  AND  V THE  CORRESPONDING  EIGENVECTOR 
RE AL*8  l I MAX, l I ,G( MAX, 1) , S21 1 1 > . VII ) , S. F ,GI » VV, S2 , SUM.E , El 
COMMON  /MATRIX/ SI  2 5, 20) ,F (20,20! ,GI I 20, 20) , VVI 100) , S2 1 1 0 , 10 ) 

CALCULATE  S»  I G 1NVERS£I*Z 

D01 I *1 , N 
001 J* l • N 
Gil  I , J ) =GI I , J ) 

CALL  D08INV(G1,N,MAX) 

0021  = 1, N 
032 J*1,N 
SC  I , J ) *0. 000 
D32K*1,N 

SII,J)*S(I,J)*GI(I,K)*2IR,J) 

IMPROVE  S THRU  ITERATION 
F-S*S-Z 

S-S-CG  INVERSE ) *F 
ND2-N 

DQIO-OI TER-1,  ND2 
D03 I =1 , N 
033  J = 1 , N 
SUM-U.ODO 
OOAK-l.N 

SUM-SUM*GC I,<)*SIK,JI 
FC I «J)*SUM-ZI I ,0) 

DOS  1 = 1 , N 
D05J-1.N 
SUM-G.ODO 
006X  = 1 , N 

SUM>SJM»GI I I,K)*F|K,J J 
SCI,J)«S|I,J)-SUM 
CONTINUE 

INITIALLY  V*(SI1,1>,  ...»  SI N»N) ) /SC  1,11,  E-l 
ITERATE:  EVeC  VV-S-V,  EVAL  El-VVCl) 

DOT  1-1 , N 

7 VID-SI  I , I ) /S  C 1 , 1 ! 

E-l. COO 

NNN»N*N 

ITER-1 

ICNT-l 

8 ITER-ITER*1 
D09I-1.N 
VVI I 1-0,000 
009J-1.N 

9 vvm-vvc  i)*sil,J)*vij) 

D010I-1.N 


10  vm-vvcti/vvm 

El-VVC 1 ) 

SJM-DABSt IE1-EI/E1) 

E-Ei 

IFISUM.GT.l. '*0-8. AND. ITER. LT.NMNIGO  TO  8 
I CNT-ICNT ♦! 

IF! ICNT.LT.5«AN0.ITER.LT «NMN)GO  TO  8 
WRITE|6,11)ITER,E,SUM 

11  F0RMATI13X,MTER«',li,'  MAX. EIGENVALUE- • ,013.6 , • ERROR-* ,013.6) 
RETURN 

ENO 
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SUBROUTINE : 
PURPOSE: 


FLOW  CHART: 


SOLVE2 

CALCULATE  MAXIMUM  EIGNEVALUE  AND  EIGENVECTOR 


CALCULATE  FIRST  (N+l)  VALUES  OF  EIGENVECTOR 
AND  MAXIMUM  EIGENVALUE 
(CALL  SOLVE!) 

i 

CALCULATE  REMAINDER  OF  EIGENVECTOR 
XLAMDA  (NP1+I)  THROUGH  XLAMDA  (NP1+NR) 


NCSL  TM-204-77 

SUBROUTINE:  SOLVE2 


DESCRIPTION:  Subroutine  SOLVE2  calculates  the  maximum  EIGENVALUE  and  EIGENVECTOR 

when  noise  is  present  on  both  input  and  output  or  the  noise  is  present 
only  on  the  output.  S0LVE2  uses  subroutine  S0LVE1  to  calculate  the 
first  (N+l)  EIGENVALUES  and  the  maximum  EIGENVALUE.  S ,'LVE2  uses 
an  interation  process  (similar  to  S0LVE1  and  S0LVE3)  to  reduce 
computation  errors. 

PROGRAM  VARIABLES: 

G GRAM  MATRIX 

GAMMA  COEFFICIENT  MATRIX 

MAX  MAXIMUM  ROWS  PERMISSIBLE 

NP1  ORDER  OF  SYSTEM+1  (N+l) 

NR  NP1-IREM 

XLAMDA  EIGENVECTOR 

Z NOISE  CORRECTION  MATRIX 


l-l 


AD-A050  206  UNIVERSITY  OF  SOUTH  FLORIDA  TAMPA  F/6  9/2 

THEORY  AND  COMPUTER  PROGRAM  FOR  TRANSFER  FUNCTION  IDENTIFICATIO--ETC (U) 
FEB  78  V K JAIN*  G J DObECK,  L J LAWDERMlLT  N61331-75-C-0012 
UNCLASSIFIED  NCSL-TM-204-78  NL 

2 of  2 I 


□ END 

DATE 

FILME 

3- 

DDC 


3-78 

DDC 
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BEST  AVAILABLE  .COPY 

, fci*  * » 


SUBROUTINE  Sr’LVE2(Z,G,GAMMA,<LAM0A«NPl»NR,M4X) 

REALMS  ZtMAX.ll.GIMAX.i) .GAMMA)  L ) , XL AMD4 ( 1 1 , SUMW , S . Z 1 , S 1 . S2 . F 
COMMON  /MATRIX /S(2  .20I.F  (25, 2"),  21  120,20 1,  SI  I10,lu»,S2tlJ,l0J 
0062 ^ 1=1 , NR 
D362C  J = 1 » NR 
ZI ( I , J )=S(NP’ ♦! ,NPl*JJ 
6 2j  S1(I,J)=ZI(I,J> 

CALL  DOR INV I Z I . NR , MAX  I 
D362il=l,NR 
0062 1J  = 1 » NP 1 
S2  ( I ,J)*'I.OOM 
0362 1K=1 , NR 

C CALCULATE  S2=G( 22)  INVERSE  * GI21) 

621  S2( I » J ) = S 2 ( I.JI+ZI ( I ,K I*GINP1 +K, J ) 

C IMPROVE  $2  THRU  ITERATION 

C F=3(22)*S2-G(21> 

C S2=S2-IG(22)  INVERSE ) *F 

P . } 

036241  TER* 1 .HR 
006221*1. NR 
00622J* 1 » NP1 
SJMw=o.r,oc 
D3623K* l . NR 

623  SUMW*SJMW+S1< I,K)*S2IK,J) 

622  FI  I ,J)=SUMW-G(NP1*I«J) 

D06  2SI =1 » NR 

00625 J = X « NP1 
SUMW=0.00O 
00626X=1 , NR 

626  SUMW=SJMW+ZI( I,K)*F(K,J) 

625  S2I I . J ) *S2 I I.JI-SUMW 

624  CONTINUE 
036271=1, NP1 
D0627J*! ,NP1 
S'JMW  = C’«  ODD 
D3628K=l,NR 

628  SJMM=SUMW+G(t.K+NPn*S2<K,JI 

627  G(  I ,J)=G( I , JJ-SUMW 
C 

C CALCULATE  XLAIMDA  1 (CALL  SOLVE  1) 

CALL  SOLVEKZ.G,  GAMMA,  XLAMDA.NPl,  MAX) 

C CALCULATE  XLAMDA  2 


036291*1, NR 
K*NP1«  I 

XLAMDA  IK)*0»'' 03 
03629J* 1 , NP1 

629  XLAMDAKI  « XL«  MD  A ( K I 2 ( I , J»*XLAM0A(  J) 
RETURN 
END 
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SUBROUTINE: 

PURPOSE: 


FLOW  CHART: 


SOLVE  3 

CALCULATE  MAXIMUM  EIGENVALUE  AND  CORRESPONDING  EIGENVECTOR 
(NOISE  ON  INPUT  ONLY) 
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SUBROUTINE:  SOLVE3 


DESCRIPTION:  SOLVE3  is  used  to  calculate  the  EIGENVECTOR  when  the  system  has 
noise  only  on  input.  S0LVE3  uses  S0LVE1  to  calculate  the  maximum 
EIGENVALUE  and  the  (N+l)  EIGENVALUES  • The  remaining  EIGENVECTOR 
elements  are  calculated  in  subroutine  S0LVE3.  S0LVE3  uses  an 
iteration  process  (similar  to  S0LVE1  an  S0LVE2)  to  reduce  computation 
errors . 


PROGRAM  VARIABLES: 


G 

GRAM  MATRIX 

GAMMA 

COEFFICIENT  MATRIX 

MAX 

MAXIMUM  ROWS  PERMISSIBLE 

NP1 

SYSTEM  ORDER  +1,(N+1) 

NR 

NP1-IREM 

XLAMBA 

EIGENVECTOR 

Z 

NOISE  CORRECTION  MATRIX 
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S JBROUT  I NE  SOLVE3IZ,G,GAMMA,XLAMOA,NP1,NR.MA<I 

RE Al*8  Z I MAX, 1) t G( MAX  < 1 ) , GAMMA! 1 » , XLAMO A ( 1 » , SUMW, S . Z I . SI . S2 ,F 

COMMON  /MATR!X/SI2.',2Ct,F(20,2C  I , Z 1 1 20, 20 J , SI  1 10, 10  I . S2 1 1C  , 10 ) , 

NPNP2*NPl*NR 

DOI 1*1 , NPl 

001 J=l,NPl 

1 ZI(1.J)>G(1.J) 

CALL  OOBINVI M »NP1 »MAX) 

D02I*1,NP1 
002 J-l , NR 
S2(  I,J»-3.0D’' 

D02X*1,NP1 

2 S2( I,J)*S2(I. JI*ZI( l,K)*GCK,J*NPl) 

D03IT£R*1,NP1 

D0AI-1.NP1 
OOAJ-l.NR 
SUMW-O* 000 
DD5K-1.NP1 

i $UMW*$JMW*GI I »K  J-S2IK  » J) 

A F(  I , JI-SUMW-GII  ,J«-NP1> 

0061*1. NP1 
006 J-l. NR 
SUMW-0.000 
D07K-1.NP1 

7 SUMW*SUMW»zm,K>*F|K,J) 

6 S2(I.J)*S2(I, J J-SUMW 

3 CONTINUE 
DOB  1*1 , NR 
009 J=1 , NR 
SUMW-O.ODO 
D09K-I.NP1 

9 SJMM*SUMri«G(T4NPl.K)*S2(K,J) 

B Gl I .J)*G( I ♦NPl , J^NPl I -SUMW 

CALL  S3LVE1(7,G,GAMMA,XLAM0A,NR,MAX) 

00101*1, NR 

10  XLAMDAINPIMI-XLAMDAI I I 
00111*1. NPl 

XLAMOAI  11*0. ‘'05 
OOt 1 J-l » NR 

11  XL A MO A (I)*XLAMDA(I)-S2(I,J) *XLAMOA( J+NP1 I 
00121*2, NPNP? 

12  XLAMOAI  n*XL»MDA(  I I /XLAMOAI  1) 

XLAMOAI  11-1. ’'DO 

RETURN  . * 

END 
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SUBROUTINE:  ZTOS 

PURPOSE:  CONVERTS  A DISCRETE  TIME  TRANSFER  FUNCTION  H(Z)  TO  A 

CONTINUOUS  TIME  TRANSFER  FUNCTION  H(S). 


FLOW  CHART: 
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SUBROUTINE  ZTOS 


DESCRIPTION:  This  subroutine  uses  the  H(Z)  transfer  function  in 

polynomial  form  and  finds  the  continuous  time  domain 
transfer  function.  Four  options  are  available  under 
the  parameter  IZTS.  If  IZTS  = 0 no  Z to  S trans- 
formation is  done.  If  IZTS  = 1 a LOGARITHMIC  trans- 
formation is  done.  When  IZTS  = 2 a PULSE  DELAYED 
transformation  is  performed.  The  fourth  option 
(IZTS  = 3)  performes  both  the  LOGARITHMIC  and  PULSE 
DELAYED  transformations. 

PROGRAM 

VARIABLES:  A NUMERATOR  POLYNOMIAL 

B DENOMINATOR  POLYNOMIAL 

DELTA  SAMPLING  INTERVAL 

N ORDER  OF  SYSTEM 

NN  ORDER  OF  NUMERATOR 

IZTS  OPTION  TO  DESIGNATE  TYPE  OF  TRANSFORMATION  DESIRED 

IZTS  - 0 PRINTS  Z DOMAIN  NUMERATOR  AND  DENOMINATOR 
AND  POLES  OF  Z DOMAIN.  (DOES  NOT  PERFORM 
Z TO  S TRANSFORMATION). 

= 1 LOGARITHMIC  TRANSFORMATION 

* 2 PULSE  DELAYED 

= 3 BOTH  LOGARITHMIC  AND  PULSE 
DELAYED  TRANSFORMATIONS  ARE 
PERFORMED. 


oonooononnnonno 
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SUBROUTINE  ZTOS 1 8# A,N, CELT*# I ITS) 
COMMON  NN 


CONVERSION  OF  A DISCRETE  TIME  SYSTEM  HU)  TO  A CONTINUOUS  TIME  SYSTEM  H(S| 


HI2)>UI1)  ♦ A!2)*ZETA  +....I/I1  ♦ 8(2)*ZETA  ♦....) 
ZETA  * L/Z 

H(SMAIl)  ♦ AI2)*S  ♦ .......  ♦ A!N*1)*$**N)/DEN0M 

DEN0M*BI1)  ♦ 8(2 )*S  ♦ ♦ 8(N*1)*S**N 

B(l)  - 1 ALWAYS 


DIMENSION  B(9)#A!9)»TEPP!20)»RRI20)#Rl(20)#CR(20)»CA!20)»CAA(20)» 
1CA1 (20) «CB<  2n) # CF( 20) ,CF1 ( 20 ) #CG( 20) 

C0MPLEX*16  CA,CAA#CA1, CB.CR. CONI, C0N2.C0NT, FAC. A1,A2.B1,B2,AA1.BBI 
XCG.CFl.CF 

REAL*8  B»A#TEMP«RR,Rl»OELTA 
CONT *0.0000 
IORP-IZTS 
NPl«N*l 
NNP1*NN*1 
A1 >0.000 
61*0.000 
DO  301-1, NP1 
IF ( t.LE.NNPl) A1*A1*A(  I ) 

30  BL>B1*B( 1 1 

WRITE! 6# 9B9)  NN.NNPl 

989  FORMAT! 10X, *NN> • , I 5,5X,*NNP1>* , IS) 

999  FORMAT!///) 

WRITE(6*999) 

WRI TE( 6# 1000) 

1000  FORMAT! • Z-DOMAIN  DENOMINATOR*) 

CALL  PRVEC(B.NPl) 

WRI TE( 6# 1001) 

1001  FORMAT! • Z-DOMAIN  NUMERATOR*) 

CALL  PRVEC(A«NP1) 

IF! IZTS.EO.O)  GO  TO  909 
IF!  IZTS.EQ.l)  GO  TO  200 
IF!  IZTS.EQ.2)  GO  TO  250 
IF!  IZTS.EQ.3)  GO  TO  200 
IF!  I ZTS.E0.4)  GO  TO  250 
200  CONTINUE 
C 

c 

C LOGARITHMIC  TRANSFORMATION 

c 

c 

C WORK  ON  NUMERATOR 

c 
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IF  I NN. EQ.O ) GO  TO  469 

CALL  POLRT I A. TEMP.NN.RR.R  I » IERI 

OOI5  I-I.NN 

15  CAI I 1-OCMPLXI RRI I1.RII11I 
00  7 I-l.NN 

7 CAI I 1- I +1.0/DELTA) *COLOGICAt It) 

IFINN.E0.N1  G0T0471 

469  CONTINUE 

DO  470  I-NNP1.NP1 
CAA III >0.000 

470  CAI I 1*0.000 

471  CONTINUE 
IFINN.EQ.OICAAI 11-1.000 

C 

C 

C NOW  THE  FIRST  NN  ENTRIES  OF  CA  CONTAIN  THE  S-DOMAIN  ZEROES  Of  NUMERATOR 
C AND  THE  REAMINING  ENTRIES  ARE  ZEROEO  OUT. 

C 

' IFINN.NE.OICALL  POLCONICA.CAA.O.N) 

C 

C 

C WORK  ON  DENOMINATOR 

c 

c 

CALL  POLRTIB.TEMP.N.RA.RI.IERl 
- 0016  I>1*N 

CRI I l-DCMPLXI  RRII).Rt(II) 

16  CFI I l-l.ODOO/CRt 1 1 

909  WRITEI6. 10021 

CALL  PRCVECICFtNI 
IFI IZTS.E0.01  GO  TO  900 
235  006  I-l.N 

6 CRI 1 1> I -1.0/DELTA }*C0L0Gt CRI I ) ) 

WRITEI6.240) 

240  FORMAT!'  LOGARITHMIC  TRANSFORMATION' ) 

WRITE! 6.9991 
WRI TEI6.2000) 

2000  FORMAT  I'  POLES  IN  S DOMAIN') 

CALL  PRCVECICR.N) 

0030001-1. N 
3000  CRI  1 1— CRI  I) 

CALL  POLCONICR.CB.O.N) 

C 

C 

C ADJUST  OC  GAIN  CONSTANT 

C 

C 

A2-CAAI 1 I 
B2-CBI1) 

FAC- IA1/B1 1*1 82/A2 1 
00  603  l-l.NNPl 
603  CAAI I 1-CAAI I) *FAC 

GO  TO  2010 
C 
C 

C OELAVEO  PULSE  INVARIANT  TRANSFORMATION 

C 

C 

C SHIFTS  NUMERATOR  COEFFICIENTS  FOR  DELAY 

C 

C 


A62 


r 
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250  CONT  = A(l) 

DO  300  I-l.N 

300  At  I)=AI !♦! » ~C0NT*B(I+1) 

AtNPU-0.0 

400  CALL  POLRT!B,TEMP,N,RR,RI,IER) 

00611-1, N 

CRt  I I-DCMPLXt  RR  1 1 ) ,R 1 1 I ) ) 

61  CFtll-l.ODOO/CRII) 

MR  I TEt  6,1002) 

1002  FORMATUX, 'THE  POLES  OF  THE  Z-DOMAIN*) 

CALL  PRCVEC tCF , N I 

C 

C 

C PARTIAL  FRACTION  EXPANSION 

C 

C 

0031-1, N 
CON  1-1,0000 
CON 2 =0.0000 
004J-1.N 

C0N2«C0N2*CR!l)*A!N-J*l> 

IFt l-J>5,4,5 

5 CONl-CONl*! l.OOOO-CRt I )*CF! J I ) 

4 CONTINUE 
3 CAC I I-C0N2/C0N1 

C 
C 

C TRANSFORMATION  OF  DENOMINATOR  ANO  NUMERATOR 

C 

C 

224  0021-1, N 

CONl-CDLOGtCRt I ) l/OELTA 

CAt I )-CAII)*CR!I )*C0N1/!CR! 1 1-1.0000) 

2 CRt  I l-CONl 

MRITEI 6,241) 

241  FORMAT t • OELAVEO  PULSE  TRANSFORMATION*) 

MRI TEt  6,999) 

226  WRITE! 6, 1004) 

1004  FORMAT! * NEGATIVE  OF  THE  POLES  IN  THE  S-OOMAIN* ) 
CALL  PRCVECtCR,N) 

MRITEI6, 1003) 

1003  FORMAT! IX, 'NUMERATOR  CONSTANTS  OF  FACTORIZED  HIS)*) 
CALL  PRCVECtCA,N) 

CALL  P0LC0N(CR,CB,0,NI 
00711-1, NP1 
71  CAA! 1 1-0.0000 
009K-1.N 

CALL  POLCON!CR,CFl,K,N) 

D09J-1.N 

9 CAA I J) -CAA! J) ♦CF1(J)*CA!K) 

CAAiNPl 1-0.0000 
2010  CONTINUE 

D0450I  = l.NPl 

CAA(I)=  CAA(I)+CONT*CB(I) 

403  MR ITEI 6, 1005) 

1005  FORMAT! • S-DOMAIN  DENOMINATOR*) 

CALL  PRCVECtCB.NPi ) 

MR  I TE( 6, 1006) 

1006  FORMAT  I*  S-OOMAIN  NUMERATOR*) 
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APPENDIX  B 

z TO  s EQUIVALENCE  TRANSFORMS 


Logarithmic  Equivalence 

Equation  (3)  on  page  4 describes  a simple  way  of  finding  an  s-domain 
transfer  function  H(s)  corresponding  to  a z-domain  transfer  function  H(z) . 

It  is  implemented  in  the  computer  program  whenever  IZTS  = 1.  The  basis 

for  this  correspondence  lies  in  the  logarithmic  mapping  s * - — £n(z)  of  the 

T 

poles  and  zeros  from  the  z-plane  to  the  s-plane.  Appropriately,  it  is  called 
the  Logarithmic  Equivalence  Transform.  The  inverse  mapping  is  similarly 
defined  and  we  note  that  the  left-hand-side  of  the  s-plane  maps  into  the 
interior  of  the  unit  circle  in  the  z-plane. 

The  primary  advantage  of  the  logarithmic  equivalence  transform  is  that 
it  preserves  the  degree  of  the  numerator  from  the  z-domain  to  the  s-domain. 
However,  it  does  not  yield  a good  degree  of  invariance  of  the  output  between 
the  continuous-time  and  discrete-time  equivalent  systems  [8].  As  a consequence, 
this  method  requires  a finer  sampling  interval  compared  to  the  'pulse  inter- 
polation' method  (described  below)  in  order  to  achieve  a satisfactory  invariance 
of  the  output. 


Leading-Edge-Pulse  Equivalence 


This  method  aims  for  invariance  of  the  output  at  the  sampling  instants. 
Strictly  speaking  this  objective  cannot  be  achieved  for  every  arbitrary  in- 
put because  the  sampled  input  signal  loses  some  of  the  information  of  the 
orginal  signal.  Suitable  restrictions  must  therefore  be  placed  on  the  class 
of  inputs  for  which  the  output  invariance  is  sought.  For  example  it  is 
assumed  that  the  bandwidth  of  the  input  signal  and  the  highest  frequency  of 
the  passband  of  the  system  are  small  compared  to  the  sampling  frequency  (say 
one-tenth  or  smaller) . Under  such  an  assumption  the  input  may  be  approximated 
by  a train  of  rectangular  pulses: 

u(t)  = u ( t ) = k | u(kA)p(t-kA) 

where  ^ 

p(t)  = for  0 < t < A 
0 

Invariance  of  the  outputs  of  a)  H(z)  excited  by  u(kA)  and  that  of  b)  H(s) 
excited  by  u(t)  can  then  be  achieved  by  equivalencing  H(z)  and  H(s)  in  the 
following  manner : 


H(z) 


n Y 
.E,  — ^ — r 

1=1  l-c^z-1 


<=> 


n 

E 

i=l 


ri 

s + Pi 


H(s)  , 


P 


i 


- i ln(a±) 


r 


i 


fjpi 

(1-c^) 


-Bl- 
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This  method  yields  a high  degree  of  invariance  between  the  outputs  of 
a)  H(z)  excited  by  u(kA)  and  b)  H(s)  excited  by  the  actual  u(t).  In  this 
respect  its  superiority  over  the  logarithmic  equivalence  method  has  been 
demonstrated  by  case  studies  on  several  Navy  vehicles.  However,  it  suffers 
from  the  disadvantage  that  the  degrees  of  the  numerators  in  the  z-domain 
and  the  s-domain  do  not,  in  general,  equal  each  other. 
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APPENDIX  C 

SOLUTION  OF  A KEY  EQUATION 


As  stated  in  Section  II  three  different  cases  for  equation  (15)  arise 
depending  upon  whether  o^,  O^,  or  both  are  nonzero.  The  solution  for  these 

different  cases  is  discussed  below. 


Case  1:  Noise  on  Both  Input  and  Output  (NSPQ  ^ 0,  NSPW  ^ 0) 

Equation  (15)  may  be  written  as 

(pi  - G_1Z) A = 0 (Cl) 

so  that  the  desired  solution  A is  the  eigenvector  of  G ^ Z corresponding 
to  its  largest  eigenvalue. 


Case  2:  Noise  on  Output  Only  (NSPQ  ^ 0,  NSPW  = 0) 


By  partitioning  the  matrix  G into  four  (n+1)  x (n+1)  blocks  and 
correspondingly  partitioning  A one  obtains 


G,.  G,„ 

I 0 

"a(1)_ 

1 

o 

1 

11  12 

-6 

II 

0 0 

a<2) 

0 

21  22 



which  is  equivalent  to  solving  the  pair 


(C2) 


[(G 


11 


“ G12G22 


LG21)  - 01]  A(1)  = 0 


(C3) 


A(2>  = 


22 


(C4) 


The  first  part  is  solved  as  a usual  eigenvalue  problem.  The  eigenvector 
A(l)  corresponding  to  the  minimum  eigenvalue  is  selected,  and,  then,  from 
the  second  equation  A (2)  is  obtained.  The  desired  parameter  vector  is 
finally  obtained  as 


A(1) 

A(2) 


(C5) 


Case  3:  Noise  on  Input  Only  (NSPQ  = 0,  NSPW  ^ 0) 

This  case  is  quite  similar  in  nature  to  case  2 above  and  is  treated 
accordingly. 

(Reverse  Page  C2  Blank) 
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